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Abstract. We have performed two-dimensional simulations of core collapse supernovae that encompass shock revival by 
neutrino heating, neutrino-driven convection, explosive nucleosynthesis, the growth of Rayleigh-Taylor instabilities, and the 
^SJ ■ propagation of newly formed metal clumps through the exploding star. A simulation of a Type II explosion in a 15 Mq blue 

^ ' supergiant progenitor is presented, that confirms our earlier Type II models and extends their validity to times as late as 5.5 

hours after core bounce. We also study a Type Ib-like explosion, by simply removing the hydrogen envelope of the progenitor 
' model. This allows for a first comparison of Type II and Type lb evolution. We present evidence that the hydrodynamics of core 

collapse supernovae beyond shock revival differs markedly from the results of simulations that have followed the Rayleigh- 
Taylor mixing starting from ad hoc energy deposition schemes to initiate the explosion. We find iron group elements to be 
synthesized in an anisotropic, dense, low-entropy shell that expands with velocities of ~ 17 000 km/s shortly after shock revival. 
, The growth of Rayleigh-Taylor instabilities at the Si/0 and (C+0)/He composition interfaces of the progenitor, seeded by the 

flow-structures resulting from neutrino-driven convection, leads to a fragmentation of this shell into metal-rich "clumps". This 
fragmentation starts already ~ 20 s after core bounce and is complete within the first few minutes of the explosion. During this 
time the clumps are slowed down by drag, and by the positive pressure gradient in the unstable layers. However, at f a; 300 s they 
^ ' decouple from the flow and start to propagate ballistically and subsonically through the He core, with the maximum velocities 

^ I of metals remaining constant at ~ 3500 - 5500 km/s. This early "clump decoupling" leads to significantly higher '*Ni velocities 

C/3 . at ; = 300 s than in one-dimensional models of the explosion, demonstrating that multi-dimensional effects which are at work 

within the first minutes, and which have been neglected in previous studies (especially in those which dealt with the mixing 
in Type II supernovae), are crucial. Despite comparably high initial maximum nickel velocities in both our Type II and our 
Type Ib-like model, we find that there are large diff^erences in the final maximum nickel velocities between both cases. In the 
"Type lb" model the maximum velocities of metals remain frozen in at ~ 3500 - 5500 km/s for / > 300 s, while in the Type II 
model they drop significantly for / > 1500 s. In the latter case, the massive hydrogen envelope of the progenitor forces the 
supernova shock to slow down strongly, leaving behind a reverse shock and a dense helium shell (or "wall") below the He/H 
interface. After penetrating into this dense material the metal-rich clumps possess supersonic speeds, before they are slowed 
down by drag forces to ~ 1200 km/s at a time of 20 000 s post-bounce. While, due to this deceleration, the maximum velocities 
of iron-group elements in SN 1987 A cannot be reproduced in case of the considered 15 Mq progenitor, the "Type lb" model 
is in fairly good agreement with observed clump velocities and the amount of mixing inferred for Type lb supernovae. Thus it 
appears promising for calculations of synthetic spectra and light curves. Furthermore, our simulations indicate that for Type lb 
explosions the pattern of clump formation in the ejecta is correlated with the structure of the convective pattern prevailing 
during the shock-revival phase. This might be used to deduce observational constraints for the dynamics during this early phase 
of the evolution, and the role of neutrino heating in initiating the explosion. 
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nucleosynthesis- 1 . Introduction 

The strong shock wave that has torn apart Sk -69° 202 and 
in consequence has given rise to Supernova SN 1987 A in 
the Large Magellanic Cloud on February 23rd, 1987, has also 
shattered the at that time widely shared belief that supernova 
explosions are essentially spherically symmetric events. The 
avalanche of observational data obtained since then from SN 
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1987 A has unambiguously demonstrated that the envelope of 
the progenitor star had substantially fragmented during the ex- 
plosion. Fast clumps of ^^Ni and its decay products had formed 
in the innermost regions of the ejecta and propagated out to the 
outer layers of the hydrogen envelope (i.e. close to the super- 
nova's photos phere) within on ly days after core collapse (com- 
pare e.g. Mitchell et al."200ll). For a bibliography of the ob- 
servational evidence see the r eviews of Hillebrandt & Hoflich 
(1989), Arnettetal. (1989a), McCrav ( 1 9931) . iNomoto et alJ 
(il994) . .Wooden (,1997.) . Muller. L1998.) . and the references 
therein. 

Spherical symmetry is not only broken in case of SN 

1987 A. This is apparently a generic feature, as is indi- 
cated by the spectra of a number of core collapse SNe that 
have been observed for more than one decade. Significant 
mixing and dumpiness of the eje cta were found fo r the 
Type II explosions SN 1995 V dFassia et all Il998l). SN 

1988 A (" SpvromiHo' 1991) and SN 1993 .T (ISDvromilidll994 
IWang&Hu .l994.: Wooslev et al.. .19941) . In case of Type lb 
supern ovae the indications for mixi ng are even older. For in- 
stance, iFilippenko & Sarg ent (1989") using [O I] observations 
have fo und that the ejecta of SN 1987 F were clumped. Even 
earlier, iHarkness et al ] (ll987l) in constructing synthetic spec- 
tra encountered the first evidence. They had to overpopulate 
excited states of He I by factors of ~ 10^ relative to local 
thermodynamic equilibrium (LTE) conditions in order to repro- 
duce the characteristic strong He lines of this supernova type . 
It was shown later by Lucv. (J99L) and iSwartz et alJ (Il993h 
that the presence of the He I lines and the implied departures 
from LTE can be understood in terms of impact excitations and 
ionizations by nonthermal electrons. The latter are presumed 
to originate from Compton-scattering of the y-rays from ""^Co 
decay and hence this hints toward significant outward mixing 
of ^^Co. Artificially mixed one-dimensional explosion mod- 
els of Type lb progenitors indeed yield good fits to observed 
Type lb spec tra and light curves dWooslev & Eastmanlll997[ 
IShi gevamaetal. 1990) . 

Such observations have stimulated theoretical and numer- 
ical work on hydrodynamic instabilities, which up to now 
has been dichotomized into two distinct classes. On the one 
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the aim to answer the question as to which extent convective 
instabilities are helpful for generating neutrino-driven explo- 
sions. Therefore, these simulations have been constrained to 
the early shock propagation phase up to about one second after 
core bounce. 

In contrast, hydrodynamic models of the late-time shock 
propagation and the associated formation of Rayleigh- 
Taylor instabilities in th e expanding mantle and envelope 
of the exploding^ star dArnett et alJ Il989h^ jPrvxell et al^l 
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the explosion mechanism. Due to considerable computational 
difficulties in resolving the relevant spatial and temporal scales 
all of these latter investigations relied on ad hoc procedures 
to initiate the explosion. This usually meant that some simple 
form of energy deposition (by a piston or "thermal bomb") 
was used to generate a shock wave in a progenitor model. 
The subsequent propagation of the blast wave was followed 
in a one-dimensional simulation until the shock had reached 
the (C-i-0)/He or the He/H interface. Then, the ID model 
was mapped to a multi-dimensional grid, and an assumed 
spectrum of seed perturbations was added to the radial velocity 
field in order to break the spherical symmetry of the problem 
and to trigger the growth of the instability. The studies of 
Yamada & Sato (1990, 1991), Nagataki et al. (19971 
and the recent 3D calculations of Hunge rford et al.l 
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who investigated the effects of metal mixing on X-ray and 
y-ray spectra formation, differ from this universally adopted 
approach by making use of parameterized, aspherical shock 
waves. 

All Rayleigh-Taylor calculations, that have been carried 
out to date for reproducing the mixing in SN 1987 A, share 
the same problem: With at most 2000 km/s their maximum 
^*Ni velocities are signific antly smaller than the o bserved val- 
ues of 3500 - 4000 km/s. IHerant & Ben3 (Il992') dubbed this 
problem the "nickel discrepancy". They also speculated that 
it should disappear when the "premixing" of the ejecta dur- 
ing the phase of neutrino -driven convect i on is taken into ac- 
cou nt. On the other hand iNagataki et alJ (Il998l) do not agree 
with IHerant & Ben3 . They claim that the nickel discrepancy is 
resolved if the supernova shock is ( mildly) aspherical. Simila r 
concl usions have bee n reached by 'Hunge rford etal] tOOA 
While iNagataki et alJ as well as ISungerfor detalJ do not rule 
out neutrino emission from the neutron star a s an explanation 
for the assumed asphe ricity of the ir shocks. iKhokhlov et all 
(Il999t) . IWheeler et all (^002) and IWang et alJ (l2002l) ques- 
tioned the neutrino-driven mechanism. Instead, they specu- 
lated on "jet-driven" explosions, which might originate from 
magneto-hydrodynamic effects in connection with a rapidly ro- 
tating neutron star 

The above controversy demonstrates the need for models 
which link observable features at very large radii to the actual 
energy source and the mechanism of the explosion. Without 
such models, interpretation of observational data requires cau- 
tion. With the present work, which is the first in a series of 
papers on non-spherical core collapse supernova evolution, 
we attempt to provide this link for the standard paradigm of 
neutrino-driven supernovae. We have performed the first one 
and two-dimensional hydrodynamic calculations of Type II and 
Type Ib-like explosions of blue supergiant stars that reach from 
20ms up to 20000 s (i.e. 5.5 hrs) after core bounce. Thus, we 
follow the evolution well beyond the time of shock eruption 
from the stellar photosphere. Our models include a detailed 
treatment of shock revival by neutrinos, the accompanying con- 
vection and nucleosynthesis and the growth of Rayleigh-Taylor 
instabilities at the composition interfaces of the progenitor star 
after shock passage. We employ high spatial resolution, which 
is mandatory to follow the episodes of clump formation, mix- 
ing, and clump propagation, by making use of adaptive mesh 
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refinement techniques. Preliminary results of this work, which 
covered the first five minut es in the evolution of a Type II su- 
pernova, were reported by Kifonidis et alJ(l2000l) . 

This paper is organized as follows: We start with an ac- 
count of the physical assumptions and numerical methods im- 
plemented in our computer codes in Sect. |2 We then describe 
our stellar model in Sect. |3l Our boundary conditions along 
with some technical details of our computational strategy are 
discussed in Sect. |4] This is then followed in Sect. |5] with our 
description of the early evolution up to 0.82 seconds after core 
bounce. To focus the discussion, we restrict ourselves to one 
basic shock-revival model, from which we start all subsequent 
calculations. To illustrate its physics as clearly as possible, we 
first present the one-dimensional case in Sect. 15. H and then the 
two-dimensional results in Sect. 15.21 In the same manner we 
discuss the subsequent evolution of a Type II supernova model 
in Sect. |6l In Sect. Q we consider the problem of a Type lb 
explosion by simply removing the hydrogen envelope of our 
progenitor. Section|8lfinally contains our conclusions. 

2. The codes 

There is as yet no numerical code that could satisfactorily han- 
dle all the computational difficulties encountered in the mod- 
elling of core collapse supernovae over the time scales that 
we are interested in in this paper Fortunately, however, such a 
code is not mandatory for studying the effects that are respon- 
sible for the formation and propagation of the nickel clumps. 
A crucial simplification arises from the fact that the phys- 
ical character of the explosion changes fundamentally from 
a neutrino-hydrodynamic to a purely hydrodynamic problem 
once the shock has been revived and launched successfully by 
neutrino heating. A simulation of the long-time evolution of 
the ejecta can therefore be split into two parts. A neutrino- 
hydrodynamic calculation that encompasses the challenges of 
modelling neutrino-driven explosions but does not require an 
extremely high dynamical range of the spatial resolution, since 
only the innermost stellar core needs to be resolved adequately. 
In contrast, in the second, hydrodynamic part, the physics of the 
problem simplifies to the solution of the hydrodynamic equa- 
tions because (to a good approximation) the central neutron star 
influences the dynamics of the ejecta only via its gravitational 
attraction and can otherwise be disregarded. It is this part of the 
calculation, however, where a high spatial resolution is essen- 
tial in order to resolve the growth and expansion of Rayleigh- 
Taylor instabilities. 

The separation into two largely independent problems 
might not be appropriate, if one is interested in an accurate de- 
termination of quantities that are influenced by the long-time 
hydrodynamic evolution of the layers near the neutron star. 
Among these quantities are the yields of r-process elements 
and the amount of fallback, which depend on the properties 
of the neutrino-driven wind as well as on the propagation of 
reverse shocks through the inner ejecta. We will not address 
these problems in the present work. Instead, we have devel- 
oped two different codes to solve each of the sub-tasks de- 
scribed above: A modified version of the hydrodynamics code 
of ljanka & Muller.Cl996,) (henceforth JM96) that includes neu- 



trino effects, and the adaptive mesh refinement (AMR) code 
A MRA, a FORTRAN im plementation of the AMR algorithm 
of lBerger & ColellaldlQSQl) . Here we discuss only the most im- 
portant features of these codes and give some details of our hy- 
dr odynamic advection sc heme. AMRA is described in length 
in lPlewa & MuIIeil(l200ll) . 

Our goal is to study the main observational consequences 
of hydrodynamic instabilities in different layers of a super- 
nova. We do not attempt to compete with the contemporary, 
highly sophisticated modelling of neutrino-driven explosions 
(for which an accurate treatment of neutrino transport is indis- 
pensable). For our purpose, a simple neutrino light-bulb algo- 
rithm, as the one of JM96, has the advantage to save enormous 
amounts o f computer time as com pared to detailled transport 
codes (see Rampp & Janka 2002, and the references therein). 
Omitting the core of the neutron star and replacing it by a 
boundary condition, that parameterizes its contraction and the 
radiated neutrino luminosities and spectra, we make use of the 
freedom to vary the neutrino properties. Thus we are able to 
control the explosion time scale and the final explosion energy 
in our simulations. By computing models with different values 
of the parameters, we can explore the physical possibilities as 
a consequence of neutrino-driven explosions. The local effects 
of neutrino heating and cooling by all relevant processes are 
treated reasonably well with our light-bulb scheme. However, 
a light-bulb description, of course, neglects all forms of back- 
reactions to the neutrino fluxes and spectra that result from 
neutrino absorption in the heating layer behind the supernova 
shock, and from neutrino emission associated with the accre- 
tion of matter onto the neutron star We have recently developed 
a gray, characteristics-based scheme for the neutrino transport 
that accounts for spectral and luminosity modifications due to 
such effects. This new treatment is currently applied in new two 
and three-dimensional simulations. The calculations in this pa- 
per, however, were still performed with the light-bulb code. 

The neutrino physics is combined with an enhanced version 
of the PPM hydrodynamics scheme (see below) and the equa- 
tion of state (EOS) of JM96. The baryonic component of this 
EOS consists of 4 nuclei (neutrons, protons, or-particles and a 
representative nucleus of the iron group) in nuclear statistical 
equilibrium (NSE). These four species are also used to com- 
pute the energy source terms resulting from nuclear transmu- 
tations. In addition to this small NSE "network", but without 
feedback to the EOS and the hydrodynamics, we also evolve 
a 14 species nuclear reaction network to approximately calcu- 
late the products of explosive nucleosynthesis, whose spatial 
distribution we wish to follow. The latter network consists of 
the 1 3 a-nuclei from "^He to ^^Ni and an additional tracer nu- 
cleus to which we channel the flow resulting from the reac- 
tion ^^Fe(a, y)^^Ni in case the electron fraction drops below 
0.490 and ^^Ni ceases to be the dom inant nucleus synthesized 
in the iron group re.g. lThielemann et al. 1996). In this way we 
can "mark" material that freezes out from NSE at conditions of 
neutron excess and distinguish it from ^^Ni whose yield would 
otherwise be overestimated. The 14 species network is solved 
for temperatures between 1 0^ K and SxIO^K. Above 8x10^ K, 
we assume that nuclei have been disintegrated to ff-particles. 
Of course, the 14 species network is a simphfication of the nu- 
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Fig. 1. Entropy (in units of ks per nucleon) in Model T3 10 at a time of 420 ms after core bounce. The white contour line encloses 
the region where the ^^Ni mass fraction is > 20%. a) Calculation performed with the original PROMETHEUS code showing 
odd-even decoupling. Note the deformation of the shock, b) Calculation performed with HERAKLES. 



cleosynthesis processes in a supernova, since it neglects impor- 
tant isotopes and production channels of the considered nuclei. 
Moreover, in our current code version the heating and com- 
position changes due to explosive silicon, oxygen, neon, and 
carbon burning as described by the 14 species network are not 
taken into account in the solution of the hydrodynamics, which 
implies an approximation of the thermodynamical history of 
the regions where nuclear burning occurs. However, the hydro- 
dynamic evolution should hardly be affected in the simulations 
presented in this paper. This can be seen by estimating the en- 
ergy release due to nuclear transmutations. In core-collapse su- 
pernovae three thermodynamic regimes can be distinguished. 
Below 2.1 X 10'^ K (the minimum temperature r equired for ex- 
plosive Ne/C burning; iThielemann et alJll996l) nuclear burn- 
ing time scales are large compared to the hydrodynamic ex- 
pansion time scale of the layers that contain the nuclear fuel, 
and no appreciable nucleosynthesis occurs. For temperatures 
2.1 X lO' K < r < 5 X 10^ K the burning time scales are com- 
parable to the expansion time scale, and the abundance and en- 
ergy changes need to be followed by solving a nuclear reaction 
network. Above 5 x lO' K NSE holds, i.e. the composition ad- 
justs immediately to temperature and density changes caused 



by the hydrodynamics. In our code, energy source terms due 
to explosive burning are not accounted for only in the second 
regime (in the NSE regime the energy source term is obtained 
from t he 4 species network). For the 15 Mq Woosley et aU 
(Il988l) progenitor that we employed (see Sect. [Sjl, we esti- 
mated that thermonuclear burning in layers with temperatures 
2.1 X 10'^ K < r < 5 X 10*^ K releases ~ 5 x 10'*'^ ergs. This en- 
ergy is small compared to the explosion energies of our models 
(about 1.8 X 10^' ergs), and is therefore negligible with respect 
to the dynamics. This may not be true, however, in case of other 
(e.g. more massive) progenitors. 

Gravity is taken into account in our neutrino- 
hydrodynamics code by solving Poisso n's equation in two 
spatial dimensions with the algorithm of iMiiller & Steinmetd 
lll995t) . To avoid the formation of transient shock oscillations 
in case general relativistic post collapse models are used as 
initial data (compare JM96), we add ID (i.e. spherical) general 
relativistic corrections to the 2D Newtonian gravitational 
potential as in Keil et al. ( 1996.) (see lRampp & Janka 2002. for 
details of the implementation and compare also to Ivan RipeJ 
Il979h . These corrections are especially important in the early 
phase of shock revival, since the deeper relativistic potential 
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well keeps the shock at significantly smaller radii than in the 
Newtonian case. However, Newtonian gravity becomes an 
excellent approximation at large distances from the neutron 
star. Therefore it can be used once the shock has emerged from 
the iron core and propagates through the oxygen shell. This 
is also the time when we map our explosion models to follow 
the late evolution with AMRA. In the AMR calculations 
Poisson's equation is solved in one spatial dimension with 
an angular average of the density. The equation of state that 
we employed for these simulations includes contributions 
from non-relativistic, non-degenerate electrons, photons, 
electron-positron pairs (the formulae for pressure and energy 
de nsity due to pairs can be derived from information given 
bv lWitti et alJll994i) and the 14 nuclei (treated as Boltzmann 
gases) that are included in the reaction network of our neutrino 
code. Separate continuity equations are solved for each of 
these nuclei in order to determine the mixing of the elements. 

The neutrino-hydrodynamics code as well as AMRA, make 
use of the HERAKLES hydrodynamics solver, an implementa- 
tion of the direct Eulerian version of the Piece wise Parabolic 
Method (PPM) of "Colella & WoodwaJ lll984 . HERAKLES 
origin ated from the PROMETHEUS code of ' Frvxell et~al] 
( Il99lh . It includes the entire functionahty of PROMETHEUS, 
but differs from its predecess or in the following a spects. The 
handling of multifluid flows of lprvxell et al.l (Il99ll) is replaced 
by the Consistent Mu ltifluid Advection scheme (CMA) of 
IPlewa & Muller'('l999V This allows for a significant reduction 
of numerical diffusion of nuclear species (to the level that is 
characteristic for the PPM advection scheme ) without suffer- 
ing f rom errors in local mass conservation (see Plewa & MiilleJ 
Il999i) . Further improvements in reducing the numerical dif- 
fusivity of the code are achieved by adopting the elaborate 
procedure for the flattening of interpolated profiles for the 
(p rimitive) state variable s that is suggested in the Appendix 
of lColella&WoodwardI £984). Additionally, the algorithm 
for the computation of the left and right input states for the 
Riemann problem is revised. The version of HERAKLES that 
is used as hydrodynamic solver in our AMR calculations has 
been further extended to achieve good computational perfor- 
mance even on small grid patches. For this purpose a new mem- 
ory interface was written, that allows for optimal pipelining 
of the one-dimensional hydrodynamic sweeps that result from 
dimensional splitting in multi-dimensional calculations, while 
due attention is paid to optimal cache reuse. With this memory 
interface the code can make efficient use of both superscalar as 
well as vector computer architectures. 

The most important modifications as compared to 
PROMETHEUS (or other standard implementations of the 
PPM scheme) are, however, the inclusion of new algorithms 
for the artificial viscosity and for r nulti-dime nsional shock de- 
tection, along with the addition of'Einfeldt'^s HLLE Riemann 
solver feinfeldt 1988 ). FoUowing Ouirk ( 1994 . 1997) we 
use the lEinfeldJ solver for zones inside strong grid-aUgned 
shocks, w hile we retain the (less diffusive) Riemann solver of 
IColella & Glaz (1985.) for all remaining grid cells. The latter 
changes are necessary to eliminate the o dd-even decouphng in- 
stability (Ouirk" 1 994 ' 1 997*; 'Liou 2000'), from which the origi- 
nal (Calslla^^-XilaZi solver suffers (.Kifonidis et al...200(]i) . This 



numerical failure shows up if a sufficiently strong shock is 
(nearly) aligned with one of the coordinate directions of the 
grid, and if, in addition, the flow is slightly perturbed. Many 
Riemann solvers allow these perturbations to grow without 
limit along the shock surface, thus triggering a strong rippling 
of the shock front and the post-shock state. In supernova sim- 
ulations, these perturbations strongly enhance the growth of 
hydrodynamic instabilities, since their amplitudes can exceed 
those of the intentionally introduced seed perturbations by sev- 
eral orders of magnitude. In case of neutrino-driven convec- 
tion this leads to faster growth of convective instabilities and 
angular wavelengths of convective bubbles which are signif- 
icantly larger than in a "clean" calculation. This is demon- 
strated in Fig. d] where we compare entropy plots of two sim- 
ulations of neutrino-driven convection that were conducted us- 
ing spherical coordinates and the neutrino parameters given in 
Sect. 15.11 The first of these (Fig. ^) was done with our orig- 
inal PROMETHEUS version and a resolution of 400 x 180 
zones and is affected by odd-even decoupling, while the second 
(Fig-ID') was obtained with the improved scheme implemented 
in HERAKLES using 400 x 192 zones. 

Odd-even decoupling and the associated carbuncle 
phenomenon (i.e. a local occurrence of odd-even decou- 
pling for shocks that are only partly aligned with the 
grid; buirklll994 1 19971 iLe Vequel Il998l) have gone largely 
unnoticed in the astrophysics literature and have marred 
almost all multi-dimensional supernova simulations con- 
ducted to date on cylindrical or spherical grids with hy- 
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3. The stellar model 



As in JM96 and iKifonidis et all ( l2000l) the initial model 
adopted for our study originated from the simulations of 
■Bruenn ( 19931 who followed core collapse and bounce o f the 
15 Mq blue supergiant progenitor of IWooslev et al.l (Il988l) . We 
have set up our shock revival calculations using his model 
WPE15 LS (180) at a time of 20 ms after core bounce. 
However, Bruenn's data set extends only to layers within 
the stellar He core and the original progenitor model of 
IWoosley~t al. ( 1988) is no longer available. To follow the late- 
time propagation of the shock in our subsequent AMR calcu- 
lation s we thu s had to reconstruct the outer envelope of the 
■Wooslev et alJ star. For this purpose we used a new 15 M© blue 
supergiant model (S. E. Woosley, private communication) that 
was smoothly joined to Bruenn's data by choosing a match- 
ing point at 8.4 X 10^ cm (M,. - 1.94Mo), i.e. within the He 
core. Table^summarizes the location of the composition inter- 
faces and the initial position of the Kg discontinuit y (that defines 
the bo undary of the iron core according to Wo osley & Weave3 
Il995h of the resulting "hybrid" 15 M© progenitor. Except for 
the He/H interface, the listed radii were taken from the post- 
collapse data of iBruena (.1993.) . 
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Since the velocity of the supernova shock and the associ- 
ated growth of Rayleigh-Taylor instabilities are very sensitive 
to the density profile, we have verified in a number of one- 
dimensional test calculations that the hybrid stellar model had 
no adverse effects on shock propagation. In these calculations 
we cut out the model's iron core and induced the explosion arti- 
ficially by depositing ~ 10^' ergs in the form of thermal energy 
at the inner boundary of the silicon shell. This procedure en- 
abled us to compare the hydrodynamic evolution of our hybrid 
model with induced explosions of newe r SN 1987 A progen- 
itor models computed by Woosl ev" et alJ ([B^97f) (for which no 
collapse calculations including neutrino transport were avail- 
able at the time we performed our simulations). No differences 
pointing to numerical artifacts were found. 



Table 1. Location of chemical interfaces in our progenitor 
model at 20 ms after core bounce. 
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Fe/Si 
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r [km] 260 
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29 800 
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M,[Mo] 1.25 


1.32 


1.50 


1.68 


4.20 



4. Computational setup 

4.1. Explosion models (hydrodynamic calculations 
with neutrinos) 

To set up a shock revival simulation we adopt spherical coordi- 
nates. We omit the innermost 0.848 M© of the core of the neu- 
tron star and replace it by the gravitational potential of a point 
mass and an inner boundary which acts as a neutrino source. 
This boundary is placed at a radius of 3.17 x 10* cm, which 
is somewhat below the electron neutrinosphere in Bruenn's 
model. The outer radial boundary is located at r = 1 .7 x 10"^ cm, 
inside the C-i-O-core of the star. Our grid consists of 400 non- 
equidistant radial zones that yield a maximum resolution of 
1.2 km at the neutron star and 200 km at the outer grid bound- 
ary. In our two-dimensional models 192 angular zones are dis- 
tributed uniformly between < < n and axial symmetry is 
assumed around the polar axis. Reflecting boundaries are im- 
posed atO -0 and 6 - n and a random initial seed perturbation 
is added to the velocity field on the entire grid with a modulus 
of 10"^ of the (radial) velocity of the post-collapse model. In 
the immediate upstream region of the accretion shock, this cor- 
responds to velocity perturbations with a modulus of ~ 5 x 10"^ 
of the local sound speed, while further outward this value de- 
creases with the decreasing infall velocities to ~ 10"^. 

In order to mimic the contraction of the cooling and delep- 
tonizing proto-neutron star, the inner boundary is moved in- 
ward during the computations, approximating the motion of 
the corresponding mass shell in Bruenn's calculations. This is 
achieved by making use of the moving grid implemented in 
our code (cf. JM96 for details). Free outflow is allowed for 
across the outer radial boundary. With this setup, the time-step 
resulting from the Courant-Friedrichs-Lewy stability condition 



Table 2. Computational setup for the 2D AMR calculation 
T310a. (see Sect. I6.2l l. The initial and final times over which 
the inner and outer grid boundaries were held fixed at ri„ and 
rout, are denoted by t; and tf, respectively. Ar is the (absolute) 
resolution and N'^^ the number of equidistant radial zones that 
would have been required to cover the entire star in order to 
achieve a resolution of Ar. Exponents are given in brackets. 



r. [s] 


ff [s] 


n„ [cm] 


'out [cm] 


Ar [km] 


Nf 


0.8(+0) 


2.3(+0) 


1.0(+8) 


4.8(+9) 


1.5(+1) 


2 599 933 


2.3(+0) 


8.0(+0) 


1.7(+8) 


1.4(+10) 


4.6(+l) 


847 804 


8.0(+0) 


3.4(+l) 


3.9(+8) 


4.3(+10) 


1.4(+2) 


280568 


3.4(+l) 


1.2(+2) 


9.9(+8) 


1.3(+11) 


4.2(+2) 


93 299 


1.2(-l-2) 


2.8(+2) 


2.3(+9) 


2.6(+ll) 


8.4(+2) 


46 650 


2.8(+2) 


6.5(-^2) 


3.9(+9) 


5.2(+ll) 


1.7(+3) 


23 282 


6.5(+2) 


1.5(-^3) 


6.8(+9) 


1.0(+12) 


3.4(+3) 


11634 


1.5(+3) 


3.6(+3) 


1.2(+10) 


2.1(+12) 


6.7(+3) 


5810 


3.6(+3) 


2.0(+4) 


2.0(-Hl0) 


3.9(+12) 


1.3(+4) 


3072 



is typically of the order of a few 10"^ s in one-dimensional cal- 
culations, and several 10 * s in two-dimensional simulations. 
Starting at 20 ms after core bounce the computations are con- 
tinued until 820 ms post-bounce, when the explosion energy 
has saturated and all nuclear reactions have frozen out. 

4.2. Models for the Rayleigh-Taylor mixing (AMR 
calculations) 

By performing test calculations we have found that, in order 
to enable a sufficiently detailed study of the growth of all 
relevant Rayleigh-Taylor instabilities, a spatial resolution of 
Ar < 10 (with R being the stellar radius) is required for 
at least some fraction of the time that it takes the shock to 
erupt from the surface of the star The above estimate holds 
for the case of a blue supergiant like Sk -69° 202. In case of 
a red supergiant an even smaller ratio Ar/R would be neces- 
sary. Without adaptive mesh refinement and/or a moving grid 
such calculations are currently unfeasible, especially in more 
than one spatial dimension. While this range of spatial scales 
appears formidable, for explicit codes like PPM an even more 
severe difficulty is posed by the vastly differing time scales in- 
volved in the problem. Even if the inner boundary is placed far 
away from the neutron star, and the gravitational attraction of 
the omitted stellar layers is taken into account by a point mass, 
the resulting Courant time step is still too small to allow one to 
follow the evolution over hours. The small time step is caused 
by fallback, which leads to high sound speeds in the central 
zones because matter that falls back toward the inner boundary 
is compressed and heated. The only remedies to this problem 
that we have found to be practical are either a coarsening of the 
grid or a shift of the inner boundary with time towards larger 
radii, where smaller temperatures are encountered. While the 
former is useful in ID calculations, where the remaining grid 
resolution that we achieved was still reasonable (see below), 
the latter proved to be the better choice for our two-dimensional 
calculations. 
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In all of our AMR simulations we have adopted spherical 
coordinates and a pre-chosen AMR grid hierarchy to which 
we mapped the data of our explosion models without adding 
any extra seed perturbations. At radii where data from the ex- 
plosion models are not available, the structural information is 
taken from the progenitor model. The inner regions of an ex- 
plosion model are cut out by placing the inner boundary of the 
AMR base level grid at a radius of 10^ cm, which is well in- 
side the neutrino-driven wind of the neutron star This baryonic 
mass flow is therefore neglected in the subsequent evolution. 
Instead we allowed for free outflow through the inner (as well 
as the outer) boundary. This is certainly a gross simplification 
for the evolution of the innermost regions of the flow during 
the first minute of the explosion (but becomes more realistic 
later when the neutrino wind is expected to cease). By ignoring 
the neutrino wind one obtains early fallback of matter through 
the inner boundary. In addition, the evolution of the reverse 
shock, which results from the slow-down of the main shock in 
the oxygen core, is affected. This reverse shock separates the 
wind from the outer ejecta at times later than ~ 500 ms after 
bounce (compare Sect. |5}. Without the pressure of the wind, 
the reverse shock propagates inward and is (partially) reflected 
at the inner boundary because the latter is not perfectly trans- 
mitting for numerical reasons. Although both the inward mo- 
tion and reflection of the reverse shock are expected to happen 
once the wind weakens (with the reverse shock being reflected 
at the steep density gradient near the neutron star "surface"), 
our approximations do not allow to model these phenomena 
reliably in space and time. Except for resulting in these limita- 
tions (which are not relevant for studying clump formation and 
propagation), the mapping from one code to the other worked 
very well and did not produce any noticeable artifacts. 

For our ID AMR calculation of Sect. 16.11 we employed 
a grid hierarchy with a base grid resolution of 512 zones and 
7 levels of refinement. We use two criteria to flag zones on a 
given level for refinement. First, we estimate the local trun- 
cation error for the (mass) density, momentum density, total 
energy density, and the partial densities of nuclei, pXi, using 
Richardson extrapolation. Whenever the error for one of these 
quantities is found to be > 10"^ (or > 10""^ in case of the pXj), 
we flag the corresponding zone for refinement. Independent of 
the truncation error estimate, we also flag zones in which den- 
sity jumps > 100% are encountered as compared to neighbour- 
ing zones. To refine between diff'erent levels, we use factors 
of 4 in resolution until level 5 is reached, where we switch to 
refinement factors of 2. With this scheme we achieve a maxi- 
mum resolution equivalent to 524288 equidistant zones. Since 
the outer boundary of the computational domain is kept at the 
stellar radius of 3.9 x 10'^ cm, this translates to an absolute 
resolution of ~ 75 km. Using this setup, we measure speed-up 
factors (as compared to using a uniform grid) that are as large 
as 10^. Still, however, the small time step resulting from the 
high sound speed in the innermost zones leads to long comput- 
ing times. We therefore remove the 7th level of the grid hier- 
archy after 1870 s of evolution and thereby reduce the maxi- 
mum resolution to 149 km. This zoning is kept until t = 4770 s, 
when also level 6 is discarded and the rest of the evolution 
followed with a maximum resolution of 298 km. It should be 




Fig. 2. Spacetime plot for the evolution of the logarithm of the 
density for our one-dimensional explosion. Model O310. The 
supernova shock is visible as the outermost discontinuity that 
extends diagonally through the plot. Around 500 ms a reverse 
shock forms in the inner ejecta that originates from the (hardly 
visible) slight deceleration of the main shock in the oxygen 
core. The reverse shock separates the ejecta from the neutrino- 
driven wind. 

noted that AMR allows one to introduce this coarsening in a 
rather straightforward fashion by manually resetting the refine- 
ment flags returned from the error estimation and zone flagging 
modules of the code. The grid can even be coarsened only lo- 
cally if required. Our ID simulation is stopped at f = 8500 s 
post bounce, 1700 seconds after the shock has emerged from 
the photosphere. 

For our 2D calculations we have adopted a diff'erent ap- 
proach, combining mesh refinement with a "zooming" algo- 
rithm for the radial coordinate. The AMR grid is chosen such 
that a base-level resolution of 48 x 12 zones and four levels of 
refinement with refinement factors of 4 in each dimension are 
employed. This results in a maximum resolution on the finest 
level corresponding to 3072 equidistant zones in radius and 
768 zones in angle. Zones are flagged for refinement by apply- 
ing the same criteria as in our one-dimensional calculation for 
each dimension. Adopting transmitting boundary conditions in 
the radial and reflecting boundaries in the angular direction, 
our computational volume is initially set up to span the domain 
(1000 < r/km < 48 000) x(O<0<k). While the grid bound- 
aries in 6 direction are kept unchanged throughout the calcu- 
lation, the zooming algorithm successively enlarges the radial 
extent of the grid according to Tabled when the shock is ap- 
proaching the outer grid boundary, rout- Whenever it is time to 
regrid, we also move the inner grid boundary, rin, away from 
the central remnant. This approach allows us to concentrate the 
computational effort in the post-shock and the Rayleigh-Taylor 
unstable regions while avoiding overly restrictive Courant time 
steps due to fallback. The times fi and ff over which the radial 
grid boundaries are kept fixed are given in columns one and two 
of Table |2] The "zooming algorithm" allows us to temporally 
achieve a radial resolution of Ar = 15 km in the unstable lay- 
ers, that is equivalent to covering the entire star with an effec- 
tive resolution of A^^*'^ - 2 599 933 equidistant zones. Even with 
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this combination of mesh refinement, grid enlargement and in- 
ner boundary movement, however, the computational load is 
still significant. The two-dimensional AMR calculation that we 
present in Sect. |6j and which follows the evolution of the mix- 
ing from 820ms to 5.5 hrs after core bounce, requires nearly 
2x10'^ floating point operations (this number has to be multi- 
plied by about a factor of 5 if no adaptive mesh is used). 

Finally, we note that the choice of spherical coordinates 
leads to numerical artifacts along the symmetry axis of the 
grid. The excessive outflows that we observe for and 
6 X n for the models in this paper are much weaker than those 
in lKifonidis et al.l (120001 Nevertheless, in order to avoid spu- 
rious effects, we have excluded a cone with an opening angle 
of 15° around the polar axis for evaluating our simulations for 
the velocity distributions of nucleosynthesis products and the 
angle-averaged extent of the mixing of chemical elements. 

5. The first second 

5.1. Evolution in one dimension 

In Table |3 we summarize the basic parameters of our one- 
dimensional explosion. Model O310, along with its two- 
dimensional counterpart. Model T310, that will be discussed in 
Sect. 15.21 As described in detail in JM96 the models are char- 
acterized by the initial values of the electron and heavy lepton 
neutrino luminosities (in units of lO^^erg/s), L|J ™d 
and the energy loss (in units of MqC^) and lepton loss parame- 
ters of the inner iron core, Ae and AFl, respectively. The latter 
parameter constrains the luminosity of electron antineutrinos 
while Ae determines the characteristic time scale for the decay 
of the luminosities that we assume to follow the simple expo- 
nential law 

L,, =L°e-'^'S(v,=Ve,Ve,v,) (1) 

where fz, is of the order of 700 ms. The neutrino spectra are 
prescribed in the same way as in JM96. 

Table 13 also lists the explosion time scale, fgxp, measured 
in ms and the explosion energy at the end of the shock-revival 
simulation, fiexp.si, in units of 10^' erg (note that the binding 
energy of the envelope has not yet been subtracted from the 
latter). We define these quantities as in JM96. The explosion 
energy is given by the sum of the gravitational, kinetic and in- 
ternal energy of all zones of the grid where this sum of energies 
is positive. The explosion time scale is defined as the time after 
the start of the simulation when the explosion energy exceeds 
lO'^'^erg. 

We adopt values for ^"d th^t give rise to en- 
ergetic explosions. Th is choice was motivated by the claim of 
lHerant&Ben3 ( ll992i) that premixing of the ^''Ni (be it artifi- 
cial or due to neutrino-driven convection) in models with high 
explosion energies of about 2 x 10^' ergs leads to a resolution 
of the "nickel discrepancy problem". In their simulations of 
Rayleigh-Taylor instabilities in the stellar envelope, they ob- 
tained final maximum ^^Ni velocities of up to 3000 km/s if they 
premixed the ^^Ni in their ID initial models out to mass co- 
ordinates of 1.5 Mq above the Fe/Si interface (i.e. thi-oughout 



75% of the metal core of the progenitor model of lNomoto et alJ 
1988). 

With the parameters listed in Table|3l shock revival by neu- 
trino heating is almost instant in our simulations. The shock 
shows neither long stagnation times nor phases of progres- 
sion and subsequent recession. Instead, it moves out of the 
iron core without noticeable delay. To illustrate matters, we 
show a spacetime plot of the evolution of the density with time 
for Model O310 in Fig. |2l from which one can infer a mean 
shock velocity of about 19 000 km/s after 0.2 s. Additionally, 
in Fig. 13 we depict the evolution of the most important hy- 
drodynamic and thermodynamic quantities as a function of ra- 
dius. The initial post-collapse profiles are plotted with thick 
lines, and the Si/O interface can be discerned by the associ- 
ated entropy step at logr [cm] = 8.8. One can easily recog- 
nize the transition of the accretion shock to an outward prop- 
agating shock with high post-shock velocities, the rarefaction 
generated by the explosion and the deleptonization and con- 
traction of the outer layers of the proto-neutron star in re- 
sponse to the contracting inner boundary. Figure |3l also shows 
the neutrino-driven wind, that is blown off the surface of the 
neutron star for times later than ~ 200 ms and that has also 
been found in the cal culations of iBurrows et alJ d 19951) and 
IJanka & MiilleJ ( Il996l) . Furthermore, a reverse shock is visi- 
ble, that forms about 500 ms after core bounce (compare Fig.|2j 
when the main shock enters the oxygen core of the star and 
propagates into layers with a somewhat flatter density gradi- 
ent, thereby decelerating slightly. Both shocks mark the bound- 
aries of a dense region that contains most of the ejecta, while 
the reverse shock separates the ejecta from the neutrino-driven 
wind. Shortly after the supernova shock has crossed the Fe/Si 
interface at 120 ms (at this time the latter has fallen in from 
its initial radius of r = 1.38 x 10** cm to 8 x 10^ cm), a high- 
density shell forms immediately behind the shock. It becomes 
more pronounced after f - 500 ms when the shock has passed 
the Si/O interface. The trajectory of the shell's inner bound- 
ary can be discerned in Fig.|2lin between the supernova and 
reverse shocks. This shell is one of the most interesting out- 
comes of our calculations. Three important physical processes 
take place in this part of the ejecta, whose mutual interaction 
has been neglected in previous attempts to model supernova ex- 
plosions. These processes are ^^Ni-synthesis, neutrino-driven 
convection, and the growth of Rayleigh-Taylor instabilities at 
the Si/O interface of the progenitor's metal core. We address 
the physical conditions for their occurrence in turn, focusing 
first on ^^Ni-synthesis. 

The left panels of Fig. |3 show the evolution of the den- 
sity and temperature distribution for Model O310 as a func- 
tion of the enclosed (baryonic) mass. The inner edge of the 
high-density shell is stationary at a mass coordinate of Mr ~ 
1.27 Mq. This mass shell marks the border between low and 
high-^e material in the expanding ejecta for times t > 120 ms 
(right lower panel of Fig.|3, i.e. for the epoch of ^^Ni forma- 
tion. Since ^''Ni synthesis requires temperatures in the range 
of ~ (5 - 7) X 10^ K and electron fractions > 0.49, we 
can expect that in Model O310 ^''Ni does not form in large 
amounts in regions with mass coordinates < 1 .27 M©, where 
Fe ^ 0.49. This is confirmed by Fig.|5j which shows the com- 
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Fig. 3. Snapshots of the density, entropy, velocity and electron fraction distribution in Model 03 10, 20 ms (heaviest Une), 220 ms, 
420 ms, 620 ms and 820 ms after core bounce. 

Table 3. Parameters of Models O310 and T310. The (baryonic) mass of the inner core, Mcoie, that has been excluded from the 
simulation, and the final neutron star mass, Mns, are given in solar masses. The initial and final radius of the inner core, R'^^^f. 
and R'^oif., respectively, are given in km, and its initial velocity of contraction, .R'^i-e, in km/s. The initial neutrino luminosities 
of different flavours, L° ^9 52' ^^'^ 52 ^^^^ ^^^^ main text), are given in units of 10^^ erg/s, and the energy loss of the 
inner core, Ae, in units of MqC^. The explosion energy, Zsexp.si, and the explosion time scale, fexp, are given in 10^' erg and in ms, 
respectively. 



Model 


Resolution 




dO 
core 


60 

core 


DCO 

core 






-^",.52 




Ae 




^exp 




O310 


400 


0.848 


31.7 


50 


18.1 


3.094 


3.684 


2.613 


0.0963 


0.0688 


1.59 


70 


1.12 


T310 


400 X 192 


0.848 


31.7 


50 


18.1 


3.094 


3.684 


2.613 


0.0963 


0.0688 


1.77 


62 


1.09 



position 820 ms after bounce, when the post-shock temperature 
has already dropped below 2.1 x 10^ K (compare Fig.|3 and 
nuclear reactions have frozen out. Large ^^Ni mass fractions 
(> 20%) are indeed only found for 1.27 Mq < M, < 1.37 Mq, 
i.e. they are confined to the dense shell. In total, 0.08 Mq of 
^^Ni are synthesized in Model O310. Of these, 0.04 Mq are 
produced by shock-induced Si-burning in the presupernova's 
silicon shell (whose inner boundary was initially located at 
1.32 Mq). The remaining 0.04 Mq stem from the recombina- 
tion of photo-disintegrated iron core matter with high electron 
fraction at mass coordinates between 1.27 Mq and 1.32 Mq. 
Between the neutron star surface and the inner boundary of the 
denseshell, i.e. for 1.14 Mq < M, < 1. 27 Mq, 0.13 Mq of mate- 
rial freezes out with Yg < 0.49. About 0.087 Mq of this matter 



end up in the neutronization tracer nucleus (denoted by X in 
Fig. O. The rest, which contains also a small contribution of 
^*'Ni, is mainly made up of a-particles, '*'*Ti, and ^^Fe, i.e. 
the products of a-rich freezeout in high-entropy (low-density) 
material. 

Table 13 gives an overview of the nuclear yields for our ex- 
plosion models (not including effects of eventual fallback). For 
a d iscussion of their acc uracy reg arding nume rical diffusion, 
see lKifonidis et alJ ( l200ll) . and Ple wa & MuUeJ (fr999-) who an- 
alyzed the effect in detail, especially for the case of '*'*Ti. With 
the spatial resolution used in the explosion models presented 
in this paper, the bulk of '^"^Ti is synthesized between about 
400 ms and 700 ms after core bounce at 1 .29 < M,. < 1 .36 Mq, 
in a layer where the abundance of a-particles decreases and 
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Fig. 4. Evolution of density, temperature, entropy and electron fraction as functions of the enclosed mass for Model 03 10 between 
t - 20 ms and t - 820 ms post-bounce. The positions of the Fe/Si and Si/O interfaces, and the Fe discontinuity (see text) are 
indicated by dotted vertical lines. 
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Fig. 5. Chemical composition of Model O310 versus mass at 
/ = 820 ms. X denotes the neutronization tracer nucleus. 



'^"Ca comes up (Fig.|5}- However, if the spatial resolution is in- 
creased, and thus numerical diffusion at the '*He/'"'Ca boundary 
is reduced, the width of the layer of '^'^Ti production via the re- 
action '*"Ca(Q', ■y)'^'^Ti is narrowed. This decreases the '^'^Ti yield 
by up to a factor of 3 (.Kifonidis et al..2001.) . 



Table 4. Final total elemental yields (in Mq) for Models O310 
and T310. Exponents are given in brackets. 





He 


'^C 


160 


-"Ne 




O310 
T310 


5.3(+0) 
5.3(+0) 


1.2(-1) 
1.2(-1) 


1.9(-1) 
1.9(-1) 


4.5(-2) 
4.5(-2) 


5.9(-3) 
6.1(-3) 










*'Ca 


44-pi 


O310 
T310 


8.6(-2) 
8.5(-2) 


2.2(-2) 
2.3(-2) 


6.2(-3) 
6.6(-3) 


6.2(-3) 
6.5(-3) 


4.6(-3) 
5.0(-3) 




^'^Cr 


52pg 




X 




O310 
T310 


3.6(-3) 
3.9(-3) 


3.7(-3) 
4.0(-3) 


8.0(-2) 
7.5(-2) 


8.7(-2) 
l.O(-l) 





The yield of Ni is much less sensitive to numerical dif- 
fusion. However, it cannot be determined to much higher ac- 
curacy than that of '^^Ti because the contribution of the region 
with Yf. < 0.49 to the ^^Ni production depends sensitively on 
the ratio of the Vg and Ve luminosities, LyJLy^, and its evolution 
with time. This ratio, which determines the electron fraction 
in regions close to the gain radius through neutrino-matter in- 
teractions ( Janka & Miiller 19 96. and references therein), must 
be regarded as uncertain in calculations like ours, that do not 
use accurate neutrino transport. With the parameters adopted 
for Model O310 we observe that Fe « 0.5 for M, < 1.27M0. 
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Fig. 7. Fractional mass of different elements contained within the velocity interval [vv, vv + dv, ] as a function of the radial velocity 
Vr in Model T310. The resolution is dv^ = 350 kms '. 




Fig. 6. Density and pressure profiles of Model O310 at t - 
820 ms. The position of the Fg discontinuity and of the Si/0 
interface are indicated by dotted vertical lines. Note the dense 
shell at r = 1.1x10^ cm, and the density and pressure gradients 
of opposite signs for 1.05 x 10^ cm < r < 1.35 x 10^ cm. 



However, given a higher Ly^ /Ly^ the opposite may occur, with 
Fe being increased in the neutrino-heated layers to values up 
to and even above 0.5, as seen in recent calculati ons including 
Boltzmann neutrino transport jjanka et alJl2003l) . If this were 
the case, the composition would be dominated by ^^Ni also for 
1.16 Mq < Mr < 1.27 M0, while only insignificant amounts 
of neutron-rich nuclei would be synthesized in these layers. 



Hence, the Ni yield might be increased by a factor of up to 
about 2 compared to Model O310. 

Just interior to the dense shell (and the ^''Ni-rich zone in 
Model 03 10), we must expect substantial convective activity to 
occur, because neutrino-matter interactions sustain a negative 
entropy gradient in the neutrino-heated material throughout the 
entire evolution up to 820 ms. As Fig. |3 shows, the entropy 
starts to rise steeply when one moves inward from 1 .27 M0 and 
reaches a maximum near the gain radius (and at times later than 
0.5 s at a position just downstream of the reverse shock). This 
negative entropy gradient will give rise to convective motions 
that, in contrast to the one-dimensional situation, will mix high- 
entropy, low-Fe, matter from the deeper layers of the iron-core 
with lower-entropy, high-Fg, material of the outermost iron core 
and the silicon shell. We discuss the implications of this in more 
detail in the next section. 

Finally, we wish to draw attention to the outer boundary of 
the dense shell. Figures|3and|6lindicate that after about 500 ms 
this boundary coincides with the star's Si/O interface. A pro- 
nounced negative density gradient forms at this interface which 
is accompanied by a positive pressure gradient (Fig.|6j. Such a 
configuration is Rayleigh-Taylor unstable even in the absence 
of gravity ( Chevalier 1976) . Its rather close spatial proximity to 
the convectively unstable layers at the inner edge of the dense 
shell is of decisive importance for the further evolution of the 
models. 
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5.2. Evolution in two dimensions 

Repeating the calculation of the one-dimensional model 03 10 
in two dimensions we obtained Model T310 that develops con- 
vective activity already 60 ms after core bounce. Fingers of 
buoyant, high-entropy material with Fg ~ 0.4, and an initial 
angular width of about 5° - 10° form quickly from the imposed 
seed perturbations in the neutrino heated layer near the gain ra- 
dius and start to rise toward the shock into the surrounding gas 
of lower entropy. While doing so they acquire a mushroom- 
capped shape and display a tendency to merge into bubbles of 
about 10° - 20° lateral width (Fig.[T). Between the bubbles, 
low-entropy material with ^ 0.498 sinks inward in tube-like 
flows. Compared to the one dimensional case, the convective 
rise of the high-entropy gas leads to a somewhat more efficient 
energy transport, and a reduced energy loss by the reemission 
of neutrinos from heated matter. Therefore Model T310 ex- 
ploded with an energy of /Sgxp - 1 .77 x 10^' erg, which is 1 1 % 
higher than in Model 03 10. The mean shock expansion veloc- 
ity is, however, almost the same as in Model O310. 

While the shock propagates through the Si-rich layers of the 
star, ''^Ni is synthesized in the high-density shell (Sect. 15. 1> . 
which is distorted at its base by the rising bubbles of high- 
entropy material. As in the one-dimensional case, freezeout 
of nuclear reactions in the deleptonized bubbles leads to high 
abundances of the neutronization tracer and to smaller yields 
for nuclei produced by a-rich freezeout. However, in between 
the bubbles, i.e. in the down-flows mentioned above, electron 
fractions are sufliciently high that ^^Ni is synthesized. Hence, 
the inner boundary of the ^^Ni shell possesses an aspherical 
shape closely tracing the inner edge of the low-entropy (high- 
density) material in the convective region, while the outer ^^Ni- 
boundary coincides with the (spherical) shock wave as long as 
the post-shock temperatures stay above 5 x lO' K. After com- 
plete silicon burning has ceased at t ^ 270 ms, an anisotropic 
nickel shell is left behind, while the spherical shock continues 
moving outward (see Fig.Q. 

The differences in the spatial nickel distribution between 
Model T310 and its ID analogue. Model O310, are rela- 
tively small, however. This is caused by the high initial neu- 
trino luminosities and their rapid, exponential decline (Eq. ^ 
adopted for our models. Both assumptions favour small explo- 
sion time scales which prevents convective bubbles to merge 
to large-scale structures, as in cases where the phase of con- 
vective overturn lasts for several turn-over times. Lowering 
the neutrino luminosities (and the explosion energies), we ob- 
tain stronger convection that largely distorts the shock wave 
by develo ping big bubbles of neutrino-heated materia l (see 
Janka&MuUer 1996 Kifonidis et al. 2000. Kifonidis et alJ 
200ll and lJanka et alJ I2OO1I for examples'). Adopting constant 
core luminosities in contrast to the exponential law of Eq. O, 
we can produce mo dels that exhib i t the v igorous boiling be- 
haviour reported bv iBurrows et alJ ill 9951) . Such cases can fi- 
nally develop global anisotropics, showing a dominance of the 
m - 0,1 = 1 mode of convection (see J anka et al- 2003; Scheck 
et al., in preparation). As a consequence, convection can lead 
to huge deviations of the spatial distribution of iron group nu- 
clei from spherical symmetry, much stronger than those visi- 



ble in Fig. n Convection can also change the nucleosynthetic 
yields of iron group nuclei, because it enhances neutronization 
by cycling high-Fg material from outer layers through the re- 
gions near the gain r adius where th e gas loses lepton number 
(Kifonidis et al. 200ll l.Tanka et al.l2 001). Due to our adopted Ve 
and Ve luminosities the gas cannot regain high values when 
it moves outward again. This leads to ^^Ni yields that are about 
6% smaller in Model T310 as compared to the ID case. Model 
03 10, while the final yield of the neutronization tracer shows 
the opposite behaviour and is 15% larger in Model T3 10 than in 
Model 03 10. These effects are comparatively small in Model 
T3 10, again mainly due to the short explosion time scale of this 
model. 

The distribution of the elements in velocity space is of sig- 
nificant importance for an understanding of the evolution of 
the supernova beyond the first second. Since Models 03 10 and 
T310 barely differ in this respect we discuss only the two- 
dimensional case. In Fig. we plot as functions of the radial 
velocity Vr and time, the fractions of the total mass of '^O, 
28 44jj ^jjjj 56]s^j jjj^j contained within the velocity in- 
tervals [vr,Vr + dv,] with dv^ - 350 kms"'. At f = 120ms, 
when the shock is travelling through the outermost layers of 
the iron core, substantial amounts of material from the silicon 
and oxygen shells are either still near hydrostatic equilibrium 
or falling towards the shock. Consequently the peak and wings 
of the mass distributions of the respective elements have zero 
and negative velocities. As the post-shock temperature begins 
to drop below 7 x 10^ K, reassembly of or-particles to heavier 
nuclei starts in the immediate post-shock region and leads to 
initially small abundances of '*'*Ti and ^^Ni, the bulk of which 
has velocities around 8 x lO-' km/s. 

Only 100 ms later (f = 220 ms) the explosion gains mo- 
mentum, and the shock propagates through the silicon shell and 
encounters substantially smaller maximum infall velocities of 
only 3000km/s. This leads to the cutoff of the high velocity 
wings of silicon and oxygen at -3 x 10^ km/s. Downstream 
of the shock, ^^Ni is synthesized as the dominant nucleus in 
regions of explosive silicon burning, while also "^Ti is built 
up in smaller amounts. The spatial neighborhood of these two 
nuclei is reflected in their velocity distributions, which both 
peak around 1.35 x 10'* km/s. However, in contrast to ^^Ni, 
some '^'^Ti also forms by a-rich freezeout in the rising convec- 
tive blobs that at this stage show the highest velocities on the 
grid. This causes the broad wing in the velocity distribution of 
"*^Ti up to 1 .7 X 10"* km/s, while maximum ^^Ni velocities reach 
1 .55 X 10"* km/s. The velocity distribution for the neutronization 
tracer shows a behaviour similar to '''^Ti. 

At 420 ms the supernova's explosion energy is still in- 
creasing and has led to maximum ^^Ni velocities as high as 
1.65 X lO'^km/s. Substantial amounts of silicon have also been 
accelerated outward. A stratification of these elements in veloc- 
ity space is beginning to emerge because the positive velocity 
gradient behind the shock (see Fig.|3} results in higher veloci- 
ties in the post-shock region as compared to gas in the deeper 
layers of the ejecta. Meanwhile, the neutrino-driven wind has 
begun to blow off material from the proto-neutron star surface 
with velocities in excess of 2.0 X lO^km/s (Fig. |3}. Traces 
of '*'*Ti in the wind, that has a composition dominated by a- 
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Fig. 8. a) Evolution of the density in our one-dimensional Type II supernova simulation Model T310 (left). The dotted lines 
indicate the positions of the Si/O, (CH-0)/He and He/H interfaces 0.82 s after core bounce, b) Shock velocity as a function of 
shock radius in IVIodel T310, and radial dependence of pr' for our progenitor model (right). Note the deceleration of the shock 
after it has crossed the (C+0)/He interface at log(r/cm) = 9.55 and the He/H interface at log(r/cm) - 10.85. The increase of the 
shock velocity for log(r/cm) > 12.2 is due to the steep drop of the density in the atmosphere of the star. 



particles and neutron-rich nuclei, can be seen at these veloc- 
ities. The distribution of '*'*Ti in velocity space is somewhat 
uncertain. We have noted in Sect. 15. II that the bulk of '*"'Ti is 
synthesized at the boundary of layers with high '*He and '*"Ca 
abundances. The location of this region coincides with the max- 
imum of the ''^Ni abundance (Fig.|5}- Thus, the peaks of both 
elements are found at comparable velocities (~ 1 .4 x 10^ km/s; 
Fig. 0. We note again, however, that the spatial distribution 
and the yield of '*'*Ti are seve rely affected by numerical dif- 
fusion (IPlewa&]VIulleilll999l) . The location of the "^Ti peak 
at f = 420 ms (Fig. might actually occur at much lower 
velocities, if numerical diffusion (and thereby '^'^Ti synthesis 
at the ^He/'^^Ca boundary) is substantially reduced. In this 
case, the contribution from the slower, deeper layers of the 
ejecta, just downstream of the reverse shock, might dominate. 
Indeed this material causes the second, smaller ^"^Ti peak near 
1.0 X 10"^ km/s at t - 820 ms (lower right panel of Fig.|7}. 

The velocities of the fastest regions that contain ^^Ni reach 
a maximum of 1.7 X 10** km/s around t = 620 ms and decrease 
slightly to 1.65 x lO'^km/s at t = 820 ms. It is obvious, that 
these ^^Ni velocities are much larger than those observed in 
SN 1987 A and that some slow-down of this material must have 
taken place later. In fact, we will show in Sect. I6.2l that exces- 
sive deceleration of this material during the subsequent evo- 
lution is the main problem for an explanation of the observed 
5*'Ni velocities in SN 1987 A. 

Summarizing our one and two-dimensional results for the 
evolution within the first second, we note that a Rayleigh- 
Taylor unstable density inversion can form at the Si/O interface 
of a blue supergiant. The quickly expanding ^^Ni is situated just 
interior to this interface and is distributed anisotropic ally in an 
inhomogeneous layer This results in a seed perturbation for the 
Rayleigh-Taylor unstable regions. We will show in the follow- 



ing sections that the magnitude of this perturbation is sufficient 
to induce significant outward mixing of ^*Ni. 

6. Beyond the first second: A Type II model 

6.1. One-dimensional evolution 

Once the shock wave has been launched by neutrino heating 
and the explosion energy has saturated, the further evolution 
of the supernova depends strongly on the density profile of the 
pr ogenitor star According to the analytic blast wave solutions 
of lSedovl ll 19591) . the shock decelerates whenever it encounters 
a density profile that falls off with a flatter slope than oc r"^ 
(i.e. for increasing pP), while for density profiles steeper than 
oc (i.e. for decreasing pr") i t accelerates. Since model s of 
supernova progenitors (see e.g. IWooslev & W eaver '1995') do 
not show a density structure that can be described by a sin- 
gle power-law, an unsteady shock propagation results that in 
turn gives rise to Rayleigh-Taylor unstable pressure and den- 
sity gradients at the composition interfaces of the star In this 
section we follow shock propagation beyond the first second 
with our AMR code in one spatial dimension to illustrate these 
effects. In order to obtain the closest possible approximation 
to the energetics of the two-dimensional case, with which we 
will compare our results in the following section, we start this 
ID calculation from an angularly averaged version of our 2D 
explosion. Model T310, (instead of the one-dimensional model 
O310). We will henceforth refer to the one-dimensional results 
of this section as Model T310. 

Figure |8^ displays the evolution of the density for this 
model. The locations of the Si/O, (C+0)/He, and He/H inter- 
faces at the start of the calculation (i.e. 0.82 s after bounce) 
are also indicated. In Fig.lSj) we show the shock velocity as a 
function of shock radius for times later than about one second 
after core bounce, along with the radial variation of pr^ for our 
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Fig. 9. Structure of Model T310, a) 20 s (left), and b) 300 s after bounce (right). Shown are (from top) the density, pressure, 
velocity, sound speed, and mass fractions as functions of radius. Note the unstable regions (shaded in grey) near the Si/O, 
(C+0)/He and He/H composition interfaces, where the density and pressure gradients have opposite signs. At f = 300 s, the 
density gradient at the inner edge of the dense shell that has formed below the He/H interface is in the process of steepening into 
a reverse shock (compare also Fig.|S^). Note also the enormous drop of the velocity in the layers that are enriched in ^^Ni, from 
6000 km/s at 20 s to 1300 km/s at 300 s. 



progenitor model. When the shock wave crosses the (C+0)/He 
interface at log(r/cm) = 9.55, 1.8 s after bounce, its propa- 
gation speed is as high as ~ 20 000 km/s. Thereafter, a rapid 
decrease of the shock velocity can be seen, since pr" increases, 
i.e. the density profile in the He core falls off with a more shal- 
low slope than oc r"^. The deceleration lasts until the blast wave 
has reached the outer layers of the He core and enters a narrow 
region of decreasing pr" around log(r/cm) = 10.7. This region 
is characterized by a rapid drop of the density just below the 
He/H interface (see Fig. |8t). Here the shock temporarily ac- 
celerates to about 9000 km/s (Fig-EJ)). Once it has passed the 
He/H interface around 80 s after core bounce, the evolution re- 
sembles that after the crossing of the (C+0)/He interface. A 
gradual deceleration of the shock to 4000 km/s occurs, until 
the blast wave enters the atmosphere of the star and finally ac- 
celerates to more than 20 000 km/s, while propagating off the 
numerical grid. 



Every time the shock decelerates, it leaves behind a positive 
pressure gradient which slows down the post-shock layers. The 
shocked material thus piles up and forms a high-density post- 
shock shell. Note that layers which are sufficiently far down- 
stream of the shock may be out of sonic contact with it, i.e. on 
the time scale of the supernova expansion these layers may not 
be reached by sound waves originating in the immediate post- 
shock region. In this case, the only way in which information 
can be mediated to the deeper layers is through a supersonic 
wave. Therefore their slow-down proceeds via a reverse shock, 
that forms at the inner boundary of the high-density shell of de- 
celerated matter and starts to propagate inwards in mass. The 
reverse shocks as well as the dense shells that form after the 
main shock has passed the Si/O, (C-i-0)/He and He/H inter- 
faces, can be clearly seen in Fig.|8t in the density structures for 
0.82 s, 20 s and 1500 s, respectively. It is also apparent that the 
density profiles for the latter two times show a striking similar- 
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Fig. 10. Logarithm of the total, time-integrated, linear growth rate in the unstable layers near the Si/O and (C-i-0)/He composition 
interfaces of Model T310 at different times. The linear growth rates increase with time and reach amplification factors > 10'^ 
within the first 300 seconds (i.e. five minutes) of the explosion. 



ity in the region that is bounded by the forward and the reverse 
shocks. 

In a one-dimensional calculation, where hydrodynamic in- 
stabilities are absent, dense shells, that have formed during ear- 
lier phases of the evolution, are preserved provided the nu- 
merical diffusivity of the employed hydrodynamic scheme is 
sufficiently small. This is demonstrated by the left panel of 
Fig.|9lwhich shows some flow quantities 20 s after core bounce. 
Note that two high-density peaks are visible between the for- 
ward and reverse shock at that time. The inner is the one that 
formed at the Si/O interface (Sect. !??!! . It contains mainly sil- 
icon and oxygen. The outer one has formed at the (C-i-0)/He 
interface and contains carbon and oxygen. In the layers of 
these shells that are shaded in gray in Fig. |5J negative den- 
sity and positive pressure gradients exist. Thus, already at this 
early time, Rayleigh-T aylor instabilities may grow in these re- 
gions ( Chevalier '1976V This is confirmed quantitatively by per- 
forming a linear stab ility analysis (see e.g. .Muller et al. 1991b : 
llwamoto et aljfl997l and the references therein; note that our 
calculations differ from these previous simulations because we 
start from a model that has a consistent explosion history, in- 
stead of having the explosion energy artificially deposited in 
the stellar core). 

In Fig. ^1 we plot for different times the total, time- 
integrated growth rate 



perturbations theory. The growth rate of the instability, cr, is 
given by 



— =exp 



[L 



crdt' 



(2) 



Pd\nPd\np 
P 



i.e. the factor by which a (small) perturbation with an initial 
amplitude would be magnified at time t according to hnear 



/5 /5 ' 
or or 

('e.g. lMulleretal]ll991bh . where we have assumed the fluid to 
be incompressible for simplicity. Within only 20 s after bounce, 
the amplification factor grows to about 10* in both unstable 
regions, which are separated by a narrow stable layer at 1 .6 Mq. 
Both growth rates are increasing until 300 seconds after bounce 
and saturate thereafter. The right panel of Fig. |5] shows that at 
t - 300 s the metal core of the star has piled up to a high- 
density zone that includes the two unstable interfaces. At this 
time, the shock has already passed the He/H interface and a 
third unstable region has formed at the outer edge of the dense 
shell that the shock has left behind at the He/H interface. The 
evolution of the growth rate in this region is shown in Fig. II II 
Alt - 300 s it is much smaller than the rates for the instabilities 
at the Si/O and (C-i-0)/He interfaces. However, while the latter 
reach their maxima at this epoch, the growth rate at the He/H 
interface still increases. It reaches about the same level as the 
growth rates at the two inner interfaces at 3000 seconds after 
core bounce. 

This behavi o ur is op posite to that found in the calculations 
of lMulleretal] l ll991bll . where the instability initially grows 
faster at the He/H interface and only after some time is sur- 
passed by the growth at the (C-HO)/He inte rface. It is very 
likely, that this result of Mul ler et alJ (Il991bl) is caused by the 
fact that they started their stability analysis as well as their 2D 
calculations 300 s after core bounce. Reliable calculations of 
the Rayleigh-Taylor growth, however, require consideration of 
the very early moments of the explosion. While the results of 
the linear stability analysis are certainly invalid quantitatively 
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Fig. 11. Logarithm of the total, time-integrated, linear growth 
rate in the unstable layers near the Si/O, (C+0)/He and 
He/H composition interfaces of Model T310 at different times. 
Comparing the amplitudes of the growth rates at the Si/O and 
(C-i-0)/He interfaces with those given in Fig. 1101 one can see 
that these have saturated after 300 s (in fact the curves for 300 s, 
1600 s, and 3000 s are not visible individually because they are 
identical within plotting accuracy). Note, however, that the rate 
for the instability at the He/H interface is still increasing be- 
tween 300 sand 3000 s. 

once non-linear growth sets in, the above conclusion is sub- 
stantiated by our t wo-dimensional calcul ations. It is also sup- 
ported by results of llwamoto et alJ ( Il997l) who found a similar 
behaviour for the growth of the instability at the (C-HO)/He and 
He/H interfaces of their 2D models for SN 1993 J. 

To our knowledge, large growth rates for Rayleigh-Taylor 
instabilities at the Si/O interface of a massive star have never 
been reported before. We can only speculate about the reason 
for the absence of this unstable layer in previous simulations. 
Either it was caused by the fact that the explosion was ini- 
tiated artificially by depositing energy in an ad hoc way, or 
it was caused by insufficient numerical resolution, or by dif- 
ferences in the structure of the progenitor stars that were em- 
ployed in these studies. The latter is obviously the case for the 
simulations of Arnett et'al] lll989bl) . Iprvxell et alJ lll99ll) and 
iMiiller et alJ ( 1199 1b) who used a stellar model at the end of 
carbon burning for their calculations. Also a combination of 
the effects listed above is possible. High-resolution studies with 
our hydrodynamic code applied to different progenitor models 
may be required to clarify this issue. We stress this point be- 
cause the two-dimensional models to be discussed in the next 
section demonstrated that the instability at the Si/O interface is 
the most crucial one for determining the subsequent evolution. 

6.2. Two-dimensional evolution 
6.2.1. Clump formation and mixing 

In contrast to the one-dimensional case, in two (and three) 
spatial dimensions the dense shells that form at the composi- 
tion interfaces can fragment into a set of smaller clumps that 



may decouple from the flow and move ballistically through the 
ejecta. Such a fragmentation actually occurs in our models due 
to the spatial proximity of the Rayleigh-Taylor unstable zones 
near the Si/O and (C+0)/He interfaces and the convective lay- 
ers. To demonstrate these effects, we present a two-dimensional 
Type 11 supernova calculation that was started from the shock- 
revival model T310. In what follows we will refer to this simu- 
lation as Model T3 10a. To illustrate its hydrodynamic evolution 
we show plots' of the mass density and the partial densities of 
the nuclei ""O, ^^Si and ^''Ni. 

Figure I12h displays the situation 4 s after core bounce, 
when the supernova shock has already crossed the (C-i-0)/He 
interface of the star The anisotropic structures that can be seen 
in this figure are strikingly similar to those visible 420 ms af- 
ter bounce (Fig.^). This indicates that the low-density (high- 
entropy) layers of the ejecta, that were heated by neutrinos, 
expand essentially self-similarly during the first few seconds 
of the explosion. However, outside of the high-entropy mush- 
rooms, and exactly as in the one-dimensional case, material 
of the Si and C-hO layers of the star piles up in two dense 
shells. Ten seconds after core bounce (Fig. 112b) also the high- 
entropy gas is affected, and the low-density mushrooms have 
been compressed to flat structures. Only 10 s later (Fig. 112b ) the 
compression has become so strong that the density distribution 
no longer resembles the flow-structures of the shock revival 
phase. A second reverse shock has formed as a consequence of 
the deceleration of the main shock in the He core (Sect.lOt. 



It is visible as the dark discontinuity at r - 1.6 x 10 cm 
in Fig. I12t . The initial mushrooms have imprinted a strong 
long-wavelength perturbation both on the dense shell behind 
the Si/O interface and at the (C+0)/He interface. Superposed 
upon this perturbation, small-scale disturbances start to grow 
along the entire Si/O interface where also about 10 cusps (not 
counting the features near the symmetry axis) begin to develop 
that are located about 20° apart. Interestingly, the positions of 
these cusps do not coincide with those of the former mush- 
rooms that are located in the regions between the cusps. The 
cusps themselves seem to be pushed by the denser material that 
constituted the former down-flows between the mushrooms. It 
is exactly in these dense regions, where ^^Ni had formed dur- 
ing the first 250 ms of the explosion. Although most of the ''^Ni 
is originally located in a layer just interior to the unstable zone 
at the Si/O interface (see Fig.ll2h and compare also Fig.|9]for 
the one-dimensional model), the dense regions with their larger 
momentum cannot be slowed down as strongly as their neigh- 
bouring matter so that they penetrate outward much farther. 
This becomes evident in the plots for t - 100 s (Fig.ll3h). The 
cusps have already grown into separate fingers that start to rise 
through the (C-i-0)/He interface and impose a long-wavelength 
perturbation onto these layers. Concurrently, the smallest scale 
perturbations that are resolved on our grid have grown to small 
mushrooms at the Si/O interface and start to mix ^^Si outwards 
while '^O is mixed inwards in between them. The growth of 
Rayleigh-Taylor fingers from the former dense Ni-"pockets" 



' Movies from this simulation are available 
http : //www . mpa-garching . mpg . de/ ~ kok/MPG 
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Fig. 12. Logarithm of the density (top) and of the partial densities (bottom) of (blue), ^**Si (green), and ^^Ni (red) in Model 
T310a at selected early epochs. Data values are coded according to the color bars given for each frame. In case of the partial 
densities, colors other than red, green and blue (resulting from the superposition of these color channels) indicate mixing of the 
composition. From left to right a) f = 4s, b) f = 10s, c) f = 20s. Note the change of the radial scale. The supernova shock is 
visible as the outermost circular discontinuity in the density plots. 




Fig. 13. Same as Fig. [HI a) f = 100 s (left), b) t - 300 s (middle), c) t - 1500 s (right). Note the change of the radial scale. 




Fig. 14. Same as Fig.IHl a) t = 5000s (left), b) t = 10000s (middle), c) t = 20000s (right). Note the change of the radial scale. 
The circular orange/black boundary in the upper right panel corresponds to the outer edge of the computational domain, which 
the supernova shock has crossed long before. 
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has important consequences for the interpretation of observa- 
tional data and will be addressed in more detail in Sect.0 

At f = 100 s the supernova shock has already passed the 
He/H interface. Furthermore, the reverse shock that had formed 
due to the slow-down of the main shock in the He core starts 
to propagate inward in radius. This reverse shock decelerates 
the innermost layers of the ejecta until it is reflected at the in- 
ner boundary and moves outward again, compressing the in- 
ner metal core of the star (Figs.lTJb and c). Five minutes after 
bounce, the metal core itself is totally shredded by the instabil- 
ity. Showing a density contrast to the ambient material of up 
to a factor of 5, the fingers have grown into the typical mush- 
room shape, which is caused by Kelvin-Helmholtz instabili- 
ties (Fig. 113b). Almost the entire metal core has thereby been 
carried through a substantial fraction of the He core, while in 
turn helium-rich gas was mixed into the metal core in "pock- 
ets" between the fingers. These low-density helium tongues can 
be discerned as the black regions around the blueish oxygen 
plumes in the abundance plots of Fig. ^] Note that there is 
a second helium-rich region closer to the center of the mod- 
els, interior to the red, green, and blue-colored zones. This re- 
gion contains material that experienced a high-entropy freeze- 
out of nuclear reactions, and thus shows high abundances of 
a-particles and nuclei like "^'^Ti. In between this ff-rich layer 
and the regions of ^^Ni dominance (the patches of dense mate- 
rial colored in red and pink in the abundance panel of Fig. 
see also Fig. |5] for the one-dimensional case) one encounters 
neutron-rich nuclei. These nuclear species participate in the 
mixing. The closer they are located to the unstable interfaces, 
the stronger their spatial distribution is affected. 

At 1500 seconds after bounce, the flow in the mixing re- 
gion has become very complex because of the interaction of 
the instabilities at the former Si/O and (C-HO)/He interfaces and 
the action of Kelvin-Helmholtz instabilities in the shear flows 
along the fingers (Fig. I13t). During this phase, and depend- 
ing on the spatial resolution, our simulations show a tendency 
for the edges of the fingers to fragment into smaller structures. 
This fragmentation appears to be more pronounced for models 
with small explosion energies. In this case, the time-scale for 
the fingers to cross the He core may increase sufficiently to al- 
low for substantial growth of fluid-dynamic instabilities at the 
boundaries of the fingers. 

While the central regions of the star are turned inside 
out, a dense shell can be seen 1500 s after bounce just be- 
low the border between helium core and hydrogen envelope 
(Fig.ll3t). It is this prominent shell that we also found in the 
one-dimensional calculation (Sect. 16. Il l as a result of the de- 
celeration of the main shock in the hydrogen envelope. Note 
that like in the one-dimensional case, the density gradient at 
the inner edge of the shell steepens into a reverse shock. In the 
one-dimensional calculation, the layers of the metal core are 
well behind this reverse shock until t - 2600 s. Only then has 
the reverse shock propagated sufliciently deep inward in m ass 
to reach the (C-HO) /He interface (see also Kane et al. and 
iMiiller et alJll99l'b) . In our two-dimensional simulation, how- 
ever, the fastest mushrooms already start to penetrate through 
the reverse shock at 1500 s after bounce (Fig. I13t). This in- 
teraction of the metal-enriched clumps with the reverse shock 



has not been pointed out in any previous study of Rayleigh- 
Taylor instabilities in Type II supernova progenitors. It has im- 
portant consequences for the velocity distribution of the ele- 
ments because it leads to a strong deceleration of the clumps. 
We will address this issue in more detail below. For the mo- 
ment we only note that after penetrating through the reverse 
shock and entering the high-density He-shell, the clumps move 
supersonically relative to the ambient medium. As a result they 
are dissipating a large fraction of their kinetic energy in bow 
shocks and strong acoustic waves. The wave fronts can be seen 
in Figs.ll4h to c which show the interaction of the clumps with 
the shell between 5000 s and 20000 s after bounce. Figure fT^ 
also suggests that during this interaction the spatial distribu- 
tion of the heavy elements inside the clumps is almost entirely 
homogenized. This causes the light blue and whitish regions 
in the abundance plots (lower panels of Figs.ll4h to c) due to 
the superposition of the single colors that were assigned to the 
different elements. The interaction with the reverse shock also 
leads to mixing of ambient helium into the clumps, with the he- 
lium mass fraction becoming comparable to that of the heavy 
elements. Before, the clumps contained only small amounts of 
helium that was admixed to them during the first ~ 300 s from 
the central zone of a-rich freezeout by the development of the 
Rayleigh-Taylor instabilities. 

To facilitate a comparison of Model T310a with one- 
dimensional work, upon which most attempts to reproduce ob- 
servations of nucleosynthetic yields and their distribution are 
based, we summarize the extent of mixing as a function of 
the enclosed mass in Fig. The left panel shows the dis- 
tribution of the mass fractions for the original presupernova 
model. At a time of 20 000 s after bounce (right panel), ele- 
ments like '^O and ^^Si, that made up the original metal core 
have been mixed almost homogeneously throughout the inner 
3.4 M0, along with the newly synthesized ^^Ni. The only nu- 
clear species that are not mixed this far out in mass are the 
neutronization tracer and, to an even lesser extent, '^'^Ti. These 
nuclei were synthesized in the innermost layers of the neutrino- 
heated ejecta very close to the collapsed core. 

Note that the mixing is confined to the former He core of 
the star (i.e. the inner 4.2 Mq). This can be understood from 
Fig. The dense He-shell, that forms just below the He/H 
interface, acts like an impenetrable wall for the metal clumps 
and shields the hydrogen envelope from becoming enriched 
with heavy elements. This result may appear somewhat sur- 
prising because, as we have shown by the linear stability anal- 
ysis (Sect. 16. it . the He/H interface of the .Woosley et al. 
star is Rayleigh-Taylor unstable. However, since we did not im- 
pose any seed perturbations in our AMR simulations and since 
the supernova shock is almost perfectly spherically symmetric 
when it crosses the He/H interface, the evolution of these lay- 
ers proceeds initially one-dimensionally. This turns out to be 
true even in models where much more vigorous neutrino-driven 
convection imposes noticeable asphericities on the shock dur- 
ing its revival. Such asphericities, however, are smoothed out 
during the subsequent propagation of the shock through the he- 
lium core. The only deviations from spherical symmetry that 
perturb the He/H interface in Model T310a are the waves that 
are excited when the metal-enriched clumps hit the inner edge 




Fig. 15. Evolution of the extent of mixing in Model T310a. Left: Initial composition. Right: Composition 20000 s after core 
bounce. The neutronization tracer is denoted with "X". Note that all heavy elements are confined to the helium core (i.e. to the 
innermost 4.2 M© of the star) even as late as 20 000 s after core bounce. 



of the dense shell below the interface at times t > 1500 s. 
However, these perturbations encounter only a moderately un- 
stable He/H interface at these late times, in accordance with 
the linear stability analysis (Fig. II 1> which shows that after 
t ~ 1500 s the integrated growth rate, cr, increases only slightly. 
The mixing at the He/H interface is therefore rather weak. The 
instability at this interface has apparently evolved into the non- 
linear regime by f » 10 000 s, showing a multitude of small and 
somewhat bent fingers. However, these do not grow apprecia- 
bly up to 20 000 s after bounce, when we stopped our 2D sim- 
ulation. At this time the flow near the interface expands nearly 
self-similarly, i.e. the small fingers move with the same veloc- 
ity as the medium between them. Note that the perturbations 
of the He/H interface that result from the accoustic waves are 
quite different from those used to initiate the instability in all 
previous studies. In the latter, seed perturbations of the order 
of 10% of the radial velocity were imposed on the entire star 
about 300 s after shock formation. With perturbations of this 
magnitude it is indeed possible to obtain strong mixing at the 
He/H interface of our progenitor. We have found such extended 
mixing in the late evolution of the hydrodynamic model dis- 
cussed in Kifonidis et al. (2000), which suffered strongly from 
odd-even decoupling, and was therefore perturbed in the radial 
velocity by ~ 20% at f = 300 s. However, even in this model 
the nickel was mixed out only to a mass coordinate of 4.5 M© 
after 10000 s. 

6.2.2. Clump propagation and velocity distribution of 
nucleosynthetic products 

In Fig. we show for times from 4 s until 20 000 s after 
bounce the logarithm of the fractional masses for several nu- 
clei that are contained within the velocity interval [v^, + dvv] 
as a function of the radial velocity. Obviously, there occurs a 
bulk deceleration of the material from velocities as large as 
~ 14 000 km/s at a time of 4 s after bounce to ~ 5000 km/s after 
50 s. This is caused by the enormous deceleration that the shock 



and the post-shock material experience after the shock has en- 
tered the He core. The deceleration and the associated com- 
pression also cause the clear tendency of the mass distributions 
of the elements to narrow down in velocity space (Fig.ll6>. 

Shortly before the shock crosses the (CH-0)/He interface 
(i.e. for times 1 s < f < 2 s) the abundance maxima of dif- 
ferent elements are well-separated in radius. Together with the 
post-shock velocity gradient (Sect. 15. 2> this results in maxima 
of the velocity distributions of different nuclei that are well- 
separated in velocity space. However, this separation disap- 
pears within the first 50 s after bounce due to the strong com- 
pression: Being squeezed into a very dense shell, different nu- 
clear species are getting very close in radius and thus they re- 
flect the velocities that are prevailing in this shell. The first nu- 
clei to become affected by the compression are '^O, ^^Si and 
^^Ni, since these are located closest to the supernova shock. A 
narrow peak forms in the velocity distribution of these elements 
already between 4 s and 10 s after bounce (Fig.ll6>. Being cen- 
tered around 8500 km/s at f = 10 s, the peak gradually recedes 
to ~ 3500 km/s at f = 50 s. The maximum velocities of ^*Ni 
and '**Si are larger than this value by about 1500 km/s, while a 
small fraction of '^O is found at somewhat larger velocities, 
still. Elements that are located closer to the inner boundary 
of the computational domain (i.e. "''*Ti and the neutronization 
tracer) show broader distributions in velocity space during the 
first ~ 10 s of the explosion, until the strong compression also 
affects their spatial distribution and leads to a homogenization 
of all profiles in Fig.ll6laround 50 s after bounce. This homog- 
enization is actually fortunate for a comparison of our results 
to observational data, since it makes the late-time velocity dis- 
tributions of nuclei in our models fairly robust against uncer- 
tainties that affect the nucleosynthesis in the neutrino-heated 
matter In particular, the velocity profiles of ""^Ni and of the neu- 
tronization tracer nucleus are finally rather similar. Hence, it is 
not very relevant for the velocity distributions of these nuclei 
that nucleosynthesis in the inner ejecta depends on unknown 
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Fig. 16. Logarithm of the fractional mass of different elements that is contained within the velocity interval [vy, Vy + dv^-] as a 
function of the radial velocity v, in Model T310a at various epochs. The resolution is dvr ~ 130km s '. Note the different scales 
for the radial velocity in the left and right panels. 



details of the neutrino luminosities and spectra, which deter- 
mine whether ^^Ni or mere neutron-rich nuclei are produced. 

Up to a time of 300 s, when the Rayleigh-Taylor insta- 
biUties are still developing, the maximum velocity of ^^Ni 
is steadily decreasing, dropping to values of ~ 3500 km/s. 
This maximum velocity is found in the mushroom caps of 
the Rayleigh-Taylor fingers. In the early stages of their exis- 
tence these "clumps" have to propagate through a rather dense 
medium. In addition they "feel" the positive pressure gradient 



in the unstable layers. Hence they are decelerated appreciably, 
but less than the expansion of the surrounding matter is slowed 
down by this pressure gradient. Around 300 s, however, the 
density contrast between the clumps and their environment has 
grown substantially. In addition, at this time the sound speed 
is about 1400 km/s in the layers of the He core through which 
the clumps propagate, while the "background" flow is expand- 
ing with ~ 3000 km/s. Comparing this to the clump velocity of 
3500 km/s we find that at f = 300 s the clumps move relative to 
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the He core with about 500 km/s, i.e. their motion relative to the 
background is subsonic and the drag they experience is rather 
small. As a result, the clumps "decouple" from the flow and 
start to move ballistically through the remaining part of the he- 
lium core, with the maximum ^^Ni velocity remaining roughly 
constant for times between 300 s and 1500 s. At approximately 
the latter time the clumps penetrate through the reverse shock 
at the inner edge of the dense shell that has formed below the 
He/H interface. After entering this shell they are strongly de- 
celerated, and the maximum velocities of all elements drop 
to ~ 1200 km/s as is visible for a time of 20000 s after core 
bounce in Fig. These values are significantly smaller than 
those observed in SN 1987 A. Note that prior to the interaction 
of the clumps with the shell, the nickel velocities of our models 
are in accordance with the velocities observed in SN 1987 A, 
i.e. a good match to the observations is prevented by the for- 
mation of the dense shell in the outer He core of our models. 
Several effects cause the strong slow-down of the clumps in the 
shell. They can be discussed by considering the expression for 
the drag force 

Fd = CD^A, (4) 

where A is the cross sectional area of a clump, p the ambi- 
ent fluid density, v the velocity of the clumps relative to the 
ambient fluid, and cq the drag coefficient. Since the clumps 
have to pass through a reverse shock they enter an environment 
with a smaller mean expansion velocity than in the layers up- 
stream of the shock. Their velocity v relative to the new back- 
ground is therefore higher than it was before they crossed the 
shock and amounts to v 1300 km/s at f = 1500 s. In addition 
p is larger in the post-shock region. Both effects increase the 
drag by increasing the ram pressure pv^. The sound speed in 
the post-shock layers is rather small (890 km/s) ait - 1500 s. 
Therefore, the clumps move supersonically through the dense 
medium, with Mach numbers M ^ 1.5. The drag coefficient 
cd of a projectile increases steeply in the subsonic/supersonic 
transition (this is the well-known "sound barrier"). The effect 
is even more pronounced if the projectile's shape is roundish 
(as in the case considered here) and not pointed. Together with 
the increase in ram pressure, the supersonic motion leads to a 
large rise of Fq, and hence to significant energy dissipation. 

The interaction of the clumps with the reverse shock re- 
sembles the well-studied problem of a (planar) shock inter- 
acting with an overdense, spherical, interstellar cloud that is 
in pre ssure equilibrium with its surroundings (see Klein et al.l 
11994 and the references therein). In this case it has been shown 
that the interaction can be divided into four stages. In the first 
phase, the blast wave strikes the cloud driving a shock into 
the cloud and a reflected shock into the (shocked) intercloud 
medium. The reflected shocks resulting from this phase corre- 
spond to the bow-shocks that we observe and which pervade 
the entire dense helium shell in our simulation between 3000 s 
and 20000 s after bounce ('Sect. l6^m i. In the second stage the 
cloud is compressed into a pancake-like structure. Indications 
of this phenomenon can be seen in Fig. I13t . The third stage 
is the reexpansion phase in which the cloud expands laterally. 
This increases its cross section. A, and thus hastens its decel- 



eration (see Figs. I14h to c). In the fourth and final phase the 
cloud is destroyed by vorticity generation due to the develop- 
ment of Rayleigh-Taylor and Kelvin-Helmholtz instabilities at 
its surface. We do not observe a complete disintegration of the 
metal clumps according to this last phase in our simulations. 
It is not clear, whether this is due to numerical reasons, due to 
the fact that we might have simply stopped our simulations too 
early for the instabilities to grow appreciably, or due to differ- 
ences in the physics of both problems. In our case the density 
contrast between a clump and the medium upstream as well as 
downst ream of the revers e shock (the former quantity is a con- 
stant in lKlein et alJl 994) is time-dependent, because the over- 
all expansion decreases the densities. Furthermore, our reverse 
shock has a Mach number M ~ I and is thus not ver y strong, in 
contrast to the shocks studied by Klein et 

ID Hill, 

for which 

M » 1 . In addition, the mushroom caps of our Rayleigh-Taylor 
fingers have typical radii rdump ~ 0.1/?s, with being the ra- 
dius of the reverse shock. Thus they d o not satisfy the small 
clump criterion (rdump ^ O.OITJJ that lKleinetalJirr994^ as- 
sumed for their study. Finally, our problem has a more com- 
plicated geometry and equation of state. Recent experiments of 
M ~ \ shocks in air interacting with a cylin drical column of 
high-density gas llFishbinell2002l: IZoldill2002h indicate that the 
gas column is strongly ablated, developing significant vorticity 
at i ts edges. However, it is not entirely destroyed. 

iKlein et all lll994h report that for the shock-cloud interac- 
tion problem, about 120 zones per cloud radius are required 
to reduce numerical viscosity in a second-order Godunov-type 
hydrodynamic scheme to the point that converged results can 
be obtained. With ~ 20 zones per clump radius our resolution is 
much lower. Resolution studies are required in order to decide 
whether our velocity distributions are numerically converged. 
The same holds for answering the question whether the clumps 
will get dispersed and mixed with the material of the "helium 
wall", or whether they can survive as individual entities. Note 
also that in the end 3D calculations will be required since the 
drag coefficient, cd, and the cross section. A, in Eq. depend 
on clump shape, which can lead to quantitative differences in 
three-dimensional as compared to our two-dimensional simu- 
lations. The Rayleigh-Taylor fingers that one finds in a two- 
dimensional calculation are comparable to genuine 3D "mush- 
rooms" only along the polar axis of the 2D grid. Along the 
equator one actually obtains toroida l structures becau se of the 
assumption of axial symmetry (see iKane et al]l20()0|. and the 
references therein). Judging from the 2D results of 'Klein e t all 
fl994) and the 3D calculations of Stone & Norman (1992) for 
the shock-cloud interaction problem, we have little doubt, how- 
ever, that independent of the dimensionality of the calculation 
the metals must slow down appreciably. We consider this to 
be a potentially serious problem for obtaining high velocities 
for the metals in the H-envelopes of Type II supernova models. 
Efficient re-acceleration mechanisms operating at later phases 
of the evolution are currently unknown. 

Actually there is no physical reason to expect an accelera- 
tion of matter (relative to the observer's frame) by the insta- 
bilities. The pressure gradient in the Rayleigh-Taylor unsta- 
ble layers of the ejecta is positive, i.e. these layers are decel- 
erated. In fact it is this deceleration that leads to the forma- 
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tion of Rayleigh-Taylor fingers: Since overdense parcels of gas 
have a larger momentum than the neighbouring material, they 
cannot be slowed down as efficiendy as the surrounding matter 
and hence they start to move outward relative to their ambient 
medium. During no time, however, the clumps are accelerated 
with respect to the observer's rest frame. The reason why we 
obtain high ^^Ni velocities in the early evolution of our models 
is the fact that the nickel clumps decouple from the (contin- 
uously decelerated) background flow and subsequently move 
ballistically through the ejecta. By decoupling from the back- 
ground, the clumps can nearly conserve their high initial ve- 
locity beyond the time of decoupling and thus move outward 
in mass until they become comoving with the outer ejecta. It 
must be noted, though, that these high initial ^^Ni velocities 
are only obtained if clump formation starts sufficiently early, 
i.e. within the first ~ 10 s after bounce, when the velocity in the 
^^Ni-rich layers is still high (compare Fig.|5^). If clump forma- 
tion is delayed by a few 100 s, the ^''Ni velocities have already 
dropped to values < 1500 km/s (Fig.|9j3) and it becomes im- 
possible to obtain high velocities by "clump decoupling" dur- 
ing the subsequent evolution. This is the reason why earlier in- 
vestigations of Rayleigh-Taylor mixing in Type II supernovae, 
which started from one-dimensional models 300 s after bounce, 
obtained only low ^^Ni velociti es. The only exceptio n were 
some of the SPH calculations of iHerant & Bend lll992h which 
showed nickel velocities of ~ 3000 km/s. As stated by the au- 
thors, this result is a direct consequence of an assumed premix- 
ing of ^*Ni throughout 75% of the metal core of the ID ex- 
plosion model of a 20 Mq progenitor, from which they started 
their 2D calculations at f = 300 s. In contrast to their approach, 
in which this (early) mixing was put in by hand, we have at- 
tempted to perform consistent and continuous high-resolution 
2D simulations of the entire shock propagation through the star 
Doing so, we find that even if one succeeds to obtain high 
metal velocities during the first minutes of the explosion, there 
is no guarantee that these velocities will remain high through- 
out the later evolution. As we discussed earlier, a main obsta- 
cle to overcome is the interaction of the metal clumps with the 
reverse shock below the He/H interface, which dramatically re- 
duces the nickel velocities in our simulations. It is not clear to 
us why this nick el deceleration is not present in the models of 

iHera nt & Be n^ 

iHerant & Bend (Il992l) reported that the energy release due 
to the radioactive decay of ^*Ni and ""^Co is inefficient to sig- 
nificantly boost the nickel velocities over the first few months. 
However, further studies of this effect, which should start from 
models like ours and use more accurate hydrod ynamic schemes 
and higher resolution than the calculations of I Herant & Bend 
lfl992) are clearly required. At present, we have to conclude 
that high final clump velocities can only be obtained if the 
clumps decouple sufficiently early from the background flow, 
and if the helium wall and the associated reverse shock do not 
form. The latter requires that the main shock does not deceler- 
ate in the hydrogen envelope, i.e. the density gradient outside 
the He core has to be steeper than oc r"^. A situation where this 
might hold is the case of a Type lb supernova explosion. The 
progenitors of Type lb supernovae are thought to be stripped 
He cores that lack a thick hydrogen envelope. In this case the 
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Fig. 17. Same as Fig. [161 but for Model T310b and times of 
300 s, 1500 s and 3000 s after bounce. 

shock directly enters the tenuous atmosphere of the star once it 
leaves the dense part of the He-rich layers. 

7. Beyond the first second: A Type lb model? 

To confirm the effect of a steeper density decline upon clump 
propagation we have removed the outer envelope of our pro- 
genitor model and replaced it by power-law profiles for the 
density and temperature of the form 



P(r) 



T(r) = Te — , 



for r > rr 



(5) 
(6) 



Herein is the radius at which the original model was cut, 
Pc and Tc are the density and temperature at that radius, and 
n = 4.25 and m = 1.32 are the adopted power-law indices for 
the density and temperature profiles, respectively. We have cho- 
sen = 5.0 X lO"* cm, i.e. we placed the cut still inside the He 
core at a mass coordinate M,- = 3.89 Mq, while the mass of the 
new envelope amounted to I.2M0. Note that this simple setup 
of a "Type lb progenitor" lacks internal consistency. It does not 
take into account that mass-loss due to winds (or binary interac- 
tion) will gradually alter the entire structure of the star during 
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Fig. 18. Logarithm of the density in Model T310b at a time of 
1500 s after bounce. 

its evolution and that the density profile at the onset of core 
collapse obtained from a consistent stellar evolutionary calcu- 
lation will therefore look different from our model. For this 
reason the calculations discussed below can only demonstrate 
effects in principle, and the simula tions should be repea ted with 
generic Type lb progenitors (e.g. IWooslev etal]ll995h to dis- 
cuss observationally relevant aspects in detail. In the context 
of the present work we refrained from doing so but wished to 
start with our explosion model T310 as described in Sect. 15.21 
and follow clump propagation in our truncated progenitor with 
AMRA until 3000 s after bounce. At this epoch the expansion 
is already self-similar, i.e. the clumps are "frozen in" and move 
with the same velocity as the ambient medium. 

The simulation in this section, to which we will henceforth 
refer to as Model T310b, was performed with a maximum res- 
olution of 1536 X 384 zones. Otherwise we used the same com- 
putational strategy as for Model T310a of Sect. 16.21 Not sur- 
prisingly, the evolution in Model T310b proceeded identically 
to the previous Type II supernova model, T310a, for the first 
50 s, i.e. the time it took the shock to reach the mass coordinate 
where we truncated the original progenitor and where the den- 
sity in the modified model now starts to decline steeply. It is 
remarkable, however, that the velocity profiles for all chemical 



elements show only minor differences between Model T310b 
(Fig. I17> and Model T310a (compare Fig. I16> even as late as 
300 s after bounce. This also indicates that the lower resolution 
used in Model T310b has only small effects on the overall dy- 
namics. However, it suppresses smale-scale fragmentation due 
to Kelvin-Helmholtz instabilities which is visible in Fig. I13t 
but not in Fig. Significant differences between both cases 
start to develop only for times later than 300 s. As we expected, 
the formation of a reverse shock and an associated dense he- 
lium shell is absent in the outer core of Model T3 10b (compare 
Fig- El to Fig. I13t). Instead of suffering from a deceleration 
in such a dense environment the metal clumps can now ex- 
pand freely into a tenuous medium, i.e. the velocity profiles 
of all chemical elements do not change any more after 300 s 
(Fig. I17>. Maximum velocities of 3500 km/s are obtained for 
^^Si, ^*Ni, '*"'Ti, and for the neutronization tracer, while '^O ex- 
pands with velocities as high as 5500 km/s. 

Figure^jdisplays the extent of the mixing (averaged over 
angle) for this model at a time of 3000 s. The newly synthe- 
sized nuclei are mixed throughout the inner 2.0 M© of the star. 
Figure Hsl should be compa red to the plots of "Model 4B" of 
IWooslev & Eastmani (Il997l) . who computed synthetic spectra 
from one-dimensional explosion models of a 2.3 Mq He core 
(originating from a 4 M© He star) and demonstrated that very 
good agreement with observed spectra of SN 1984 L near max- 
imum light could be achieved. This was only the case, how- 
ever, if they artificially mixed ^^Ni into the helium rich lay- 
ers assuming an exponential decline of the ^^Ni mass frac- 
tion, X(^"'Ni), with the enclosed mass. The value of X(^''Ni) 
decreased below 10"^ only a t a mass coordinate of ~ 2 M(7i. 
Different from Model 4B of IWooslev & Eastma ill 9971) our 
Model T310b (Fig. I19> shows a plateau-like distribution with 
high nickel mass fractions (logX(^^Ni) ^ -1.0) up to the 
outer edge of t he nickel-enriched co re. The latter coincides 
with that of the IWooslev & Eastmani model. Whether this is 
sufficient to cause an adequately large flux of ionizing y- 
photons in the He envelope in case of our more ma ssive star 
(5.1 Mq instead of 2.3 Mq for IWooslev & Eastmaiil 's "Model 
4B") can only be decided by calculating detailed spectra. Our 
purpose here is only to demonstrate that multi-dimensional hy- 
drodynamic models, with a consistent treatment of the early 
phases of the explosion, can yield strong mixing of the metal 
and He core that is a prerequisite to produce a good match 
between calculated an d observed spectra and hght curves 
of Type lb supernovae (Shigevama et al. 'l990'; 'Hachisu e t alJ 
,lg94: Woos lev & Eastman 1937) . Hachisu et al. (.1994) have 
shown that in order to obtain a sufficient amount of mix- 
ing in multi-dimensional simulations with artificially triggered 
Rayleigh-Taylor instabilities, the amplitude of the seed pertur- 
bations must exceed 5% of the radial expansion velocity at the 
time the shock crosses the unstable interfaces in the metal core. 
Our models demonstrate that such a degree of perturbation can 
be naturally obtained as a result of neutrino-driven convection 
during the first second of the explosion. 

A most remarkable feature of Model T310b is the fact that 
the number of Rayleigh-Taylor fingers at times t > 1500 s is 
correlated with (indeed it is essentially equal to) the number of 
down-flows that developed in the convective layer of neutrino 




heating behind the shock (compare Fig. ^] to Fig. QJi). This 
correlation, although initially present, is destroyed in Model 
T310a, where the interaction of the clumps with the reverse 
shock and the dense helium shell results in a stronger non- 
linear evolution with substantial vorticity generation and mix- 
ing. As the structure of down-flows and neutrino-heated bub- 
bles characterizes the mode of the convective pattern that pre- 
vails within the shock-revival phase, it carries important in- 
formation about the mechanism that initiates the explosion. It 
would be interesting if our finding is generic, and if it also 
holds for cases where the pattern of neutrino-driven convec- 
tion develops much larger modes than those obtained in our 
simulations. At present, however, one must be cautious not to 
overinterpret our result, which is based on a single 2D model. 
A confirmation of the existence of a tight correlation between 
the patterns of the early post-shock convection and the final in- 
stabilities in the mantle of the exploding star will probably re- 
quire well-resolved 3D calculations and parameter studies with 
difi^erent explosion and progenitor models. 

8. Conclusions 

We have presented a study of hydrodynamic instabilities in a 
Type II and a Type Ib-like supernova model that has improved 
upon earlier simulations by starting the explosion by means of 
neutrino heating behind the stalled supernova shock and by 
following the hydrodynamics of this explosion well beyond 
shock-breakout from the stellar photosphere. The consistency 
of the hydrodynamic evolution that was thereby achieved al- 
lows for the following conclusions. 

- We find newly synthesized ^^Ni and the products of in- 
complete silicon and oxygen burning to be initially located 
in a high-density (low entropy) inhomogeneous shell that 
forms behind the outward sweeping shock within the first 
~ 500 ms of the evolution. Depending on the explosion 
time scale and the strength of the neutrino-driven convec- 
tion this shell is either markedly or weakly deformed by 
rising, high-entropy bubbles of neutrino-heated matter. The 



angular extent of these bubbles is larger for long explosion 
time scales (usually leading to smaller explosion energies 
of our progenitor) and large seed perturbations. Besides be- 
ing formed by Si-burning ^''Ni is also produced in signifi- 
cant amounts in the neutrino-heated high-entropy bubbles 
if their neutronization is moderate. The latter is sensitive to 
the Ve and Vg luminosities and spectra imposed at the inner 
boundary in our simulations. 

- The outer boundary of the ^''Ni-rich shell is located close 
to the Si/0 composition interface where a negative density 
and a positive pressure gradient exist. This state is unstable 
to the Rayleigh-Taylor instability. 

- For an explosion energy around 1 .8 x 10^' ergs, ^''Ni is born 
with velocities < 17 000km/s. These velocities are much 
larger than those typically observed for the metal clumps 
in supernovae, and indicate that a substantial deceleration 
of this material must occur when it penetrates through the 
overlying layers of the star. 

- Three Rayleigh-Taylor unstable regions dev elop in the ex- 

flosion of the 15 Mq progenitor model of I Wooslev et alJ 
19881) : one at the Si/O interface, a second one at the 
(C-i-0)/He interface and a third one at the He/H interface 
of the star 

- Seeded by the perturbations induced by neutrino-driven 
convection, Rayleigh-Taylor mixing at the Si/O interface 
sets in already about 20 seconds after core bounce and leads 
to a fragmentation of the ''''Ni-enriched shell with the for- 
mation of fingers that distort also the (CH-0)/He interface 
farther out. 

- Within onl y about five minutes af ter bounce the entire metal 
core of the lWooslev et alJ ill 988 ') star is shredded. 

- Dense "bullets" and clumps of ^''Ni, ^^Si and ""O-rich ma- 
terial decouple from the flow and start to propagate ballis- 
tically through the stellar helium core. 

- During this phase, metal velocities as high as 3500- 
5500 km/s are observed, and the pattern of the Rayleigh- 
Taylor structures carries information about the geometry 
of the neutrino-driven convection from the shock-revival 
phase. 
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The effects discussed above could not be seen in previous 
studies of Rayleigh-Taylor instabilities in core collapse super- 
novae where the explosion was started by ad hoc energy de- 
position (using pistons or thermal bombs). Especially the high 
^^Ni velocities, the strong mixing triggered at the Si/O inter- 
face, and the very short growth time-scale of the instabilities 
at the Si/O and (CH-0)/He interfaces (minutes as compared to 
hours) are in clear contrast to what has been reported in pre- 
vious studies, especially in those which dealt with the mixing 
in Type II supern ovae. Of all e arlier works, only the Type lb 
st udies of Hachis u et alJ ( 1 1994 and the SN 1993 J simulations 
of llwamoto et alJ ( 1 1997* ) agree qualitatively with the evolution 
that we observe for the first minutes of the explosion. Both of 
these studies found, as we do, that Rayleigh-Taylor instabilities 
at the (C-i-0)/He interface of the progenitor stars grow within 
minutes after core collapse. However, they did not observe the 
instability at the Si/O interface. Moreover, they had to assume 
large seed perturbations to induce the mixing, whereas in our 
case the perturbations are a natural consequence of neutrino- 
driven convection. 

The crucial difference between our Type II model T310a 
and all other Type II supernova simulations (except for the 
work of Iwamoto et al. 1997) is the fact that our metal core 
starts to fragment within only about 20 s after bounce, when the 
velocity in the ^^Ni-rich layers is still high (see Fig.|9^). This 
leads to an early decoupling of the clumps from the background 
flow, enabling the clumps to conserve the high velocity that 
the flow possesses at early times. Not allowing for early clump 
formation and decoupling, as it has been done in nearly all pre- 
vious multi-dimensional simulations of Type II supernovae by 
starting the calculations hundreds of seconds after core bounce 
from ID initial models, has the effect that the velocities of the 
^^Ni-rich layers have dropped to values < 1500 km/s at the start 
of the multi-dimensional calculation (Fig.|5J))- It therefore be- 
comes impossible to obtain high velocities by the Rayleigh- 
Taylor instabilities during the subsequent evolution. 

In fact, the situation becomes even worse at later times. We 
find that our Model TBlOa is initially successful in yielding 
high ^*Ni velocities: Once the ^^Ni clumps in this model have 
decoupled from the flow, they propagate ballistically and sub- 
sonically through the helium core. Due to their subsonic mo- 
tion, they do not dissipate much of their kinetic energy and the 
maximum ^^Ni velocities remain high at 3500 km/s. However, 
this changes when the clumps reach the outer helium core and 
encounter the reverse shock and the dense (Rayleigh-Taylor un- 
stable) "wall" that are left behind by the main shock below the 
He/H interface. When the clumps penetrate through the reverse 
shock into this wall, 1500 s after core bounce, they are decel- 
erated due to supersonic drag to velocities < 1500 km/s. An 
interaction of the clumps with the reverse shock i s also present 
in data o f previous Type II supernova simulations jMiiller et alJ 
Imib'a"). However, its importance for clump propagation was 
not recognized because 

- due to the late phase when previous multi-dimensional cal- 
culations were started, clump formation at the (C-i-0)/He 
interface was only about to set in at f » 2000 s, when the re- 



verse shock had propagated sufficiently far inward in mass 
to reach this interface, 

- the initial clump velocities and the densities in the shell 
downstream of the reverse shock were low, leading to less 
pronounced effects during the interaction, 

- and the employed numerical resolution was insufficient to 
resolve the non-linear waves that characterize this interac- 
tion. 

A final difference between Model T310a and previous 
Type II supernova simulations is the fact that we do not ob- 
tain strong mixing at the He/H interface. Unless we add artifi- 
cial perturbations to these layers, as it was done in all previous 
studies, the mixing remains very weak. This might indicate that 
neutrino-driven convection is not able to provide perturbations 
that are sufficiently strong to cause the large-scale mixing that 
was observed in SN 1987 A. However, final conclusions regard- 
ing this point are not possible before Rayleigh-Taylor simula- 
tions have been performed which start from a globally aspher- 
ical situatio n, as seen in recen t calculations of neutrino-driven 
explosions Jjanka et alJ20 03^ Scheck et al., in preparation). 

Taken together, the effects revealed by our simulations 
demonstrate that the evolution proceeds multi-dimensionally 
from the earliest moments. Convective instabilities during the 
shock-revival phase cannot be neglected. They lead to very 
early growth of Rayleigh-Taylor instabilities and thus to high 
initial clump velocities what in turn determines the character 
of the late-time hydrodynamic evolution. By starting from ID 
initial models hundreds of seconds after core bounce, all earlier 
investigations of Rayleigh-Taylor mixing in Type II supernovae 
have therefore missed essential parts of the physics of the prob- 
lem. 

Supernova 1987 A still remains an enigma, though. Given 
the (late) deceleration of the clumps below the He/H inter- 
face in Model T310a, we wonder whether the conflicting max- 
imum ^^Ni velocities between our models and the observa- 
tions point to deficiencies of the "standard" blue supergiant 
progenitors, or to missing physics in our shock-revival simu- 
lations. There are indications that the progenitor of SN 1987 A 
might have been the result of a merger of two smaller stars 
jPodsiadlowskil 1 99^ and that therefore its structure might not 
be accounted for reliably by "standard" stellar evolutionary 
calculations. Nume rical studies of such mergers have recently 
been performed bv llvanova et alJ ll2002l) . Although their mod- 
els show density profiles in the hydrogen envelope that look 
similar to the one of our progenitor (N. Ivanova, private com- 
munication), it remains to be seen in future simulations whether 
for these models the evolution differs from what we have re- 
ported here. 

Sources of uncertainty in our calculations that may af- 
fect the velocities of chemical elements are the still lim- 
ited numerical resolution and especially the differences be- 
tween two-dimensional versus three-dimensional hydrodynam- 
ics. The smaller drag that clumps experience in three dimen- 
sions ( Kane et al. 2000) may lead to higher initial clump ve- 
locities after the fragmentation of the metal core and thus to a 
farther penetration of the clumps through the helium "wall". 
This might result in higher final ^^Ni velocities than in two 
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dimensions. An extension of our calculations to three dimen- 
sions with a similar resolution as in the 2D case appears hardly 
feasible in view of present computer resources. Fur ther insight 
might , however, be gained by la ser experiments (Robev et al.l 
l200lt iKane etaljEooH brake e t al. 2002; Klein et al. 2003, 
and the references therein), and hydrodynamic code validation 
experi ments as those described by .Fishbine, (i2002) and Zoldi 
ll2002h before well-resolved 3D hydrodynamic simulations will 
become available. 

There has been much speculation about large-scale 
anisotropies caused by jets in the explosion of SN 1987 A 
jKhokhlov et al." 1999': ' Wheeler e t al."2002; 'W ang et alJ20 0?). 
However, convincing and consistent (MHD) calculations of the 
formation and propagation of such jets and an associated ex- 
plosion have not yet been performed. It is unclear whether in 
this scenario nickel velocities can be obtained that are in agree- 
ment with those observed in SN 1987 A. An initial anisotropy 
of the explosion, that might lead to high initial nickel velocities, 
would have to work against the same effects in the envelope that 
we find to be important in our simulations. The interaction with 
the reverse shock and the dense helium shell might not only 
slow down fast metal clumps, it might also reduce an initial 
anisotropy of the ejecta. Judging from the experience gained 
from our calculations, we doubt that the spatial resolution in 
simulat ions of anisotropic explosions that have been perfo rmed 
to date dNagataki et alll"997rfl998tE^khlov et al.ll 999h was 
sufficient to study these effects reliably. The origin o f the pro- 
late deformation of the ejecta of SN 1987 A dWang et al._2002 ) 
must therefore be regarded to be unknown. In fact, there is cur- 
rently neither observational nor theoretical evidence that un- 
ambiguously demands the assumption of a "jet-driven" explo- 
sion. The deformation could as well be linked to an anisotropic 
initiation of the explosion due to neutrino heating (' janka et all 
I2OO3; Scheck et al., in preparation). On the other hand it might 
also be caused by stellar rotation as a result of a merger history 
of the progenitor 

While clearly substantial work is requked in case of 
SN 1987 A, our "Type lb" Model T310b is in reasonably 
good agreement with observations of extragalactic Type lb 
supernovae. Due to the absence of dense shell and reverse 
shock formation in the He core of this model, the metal-rich 
clumps are not decelerated once they start to propagate bal- 
listically through the ambient gas and the final metal veloc- 
ities of 3500 - 5500 km/s are sufficiently high. In addition, 
the extent of the mixing in this model is comparable to what 
[Woosley & Eastman ( 1997) had to assume in their "Model 4B" 
to model the spectrum of SN 1984L. This result is very en- 
couraging. It indicates that the interaction of neutrino-driven 
convection with the Rayleigh-Taylor instabilities in the stellar 
mantle is able to account for important aspects of the mix- 
ing of different elements in Type lb supernovae which are 
known from spectral and light curve calculations for more than 
a decade. 

Another interesting aspect of Model T310b is the fact that 
without the strong non-linear interaction of the clumps with the 
reverse shock, perturbations originating from the shock-revival 
phase are reflected in the final flow structures of the Rayleigh- 
Taylor instabilities. This may enable one to deduce information 



about the pattern of the neutrino-driven convection from obser- 
vations of the distribution of metals in the ejecta of Type lb 
supernovae. The possibility to obtain such a unique piece of 
evidence from ejecta properties is exciting. However, before 
such conclusions can be drawn, our results need to be ascer- 
tained by making use of a larger number of shock-revival simu- 
lations with different properties and a greater range of (Type lb 
and Type Ic) progenitors. We believe that in the end very well- 
resolved 3D simulations will be required to prove the existence 
of such a correlation. If confirmed it might provide us with a 
means to probe the physics linked to the explosion mechanism 
even for extragalactic supernovae. 

Finally, we emphasize that clump deceleration due to the 
interaction with the reverse shock might be crucial for a cor- 
rect interpretation of observations of the different supernova 
types. It can cause different final metal clump velocities even 
if the initiation of the explosion proceeds similarly in differ- 
ent stars. The latter must actually be expected for the neutrino- 
heating mechanism, because the post bounce models of super- 
nova cores are rather similar, regardless of whether they orig- 
inate from Type II or Type lb progenitors (Rampp et al., in 
preparation). On the other hand, clump deceleration depends 
on the structure of the envelope of the progenitor and a se- 
quence is conceivable where its importance varies with the 
type of the supernova, being strongest in (some fraction of?) 
Type II events and weaker or absent in Type lb and Ic su- 
pernovae. Unfortunately, observational data of metal clump 
velocities are currently sparse and hence they do not allow 
one to test this hypot hesis. Photosph eric velocities, as those 
published recently by 'Hamu^ ll2003h . are not very helpful in 
this respect because they only probe the outer, faster expand- 
ing layers of the ejecta. To enlarge the data base recourse 
must be made to measurements of velocities in supernova rem- 
nants. Clumps movin g with up to ~ 6000k n i/s were, however, 
observed in Cas A jvan den Berghl fl97 1[ iThorstensen et all 
'2001') , which is probably th e remnant of a Type lb explo- 
sion jFesen&Rec keH ll99ll) . and clump velocities of up to 
~ 6000 km/s and ~ 4000 km/s were measured shortly after 
the explosion in the Type lb SN 1987 F ("Filippenko & SargentI 
■ 1989) and the Type lib SN 1993 J ( SpvromiUo 1994.) . respec- 
tively. All of these objects are assumed to be connected to 
progenitors with small or lacking hydrogen envelopes. On the 
other hand .Aschenbach 12002) and Aschenbach et all (Il995l) 
deduce low mean expansion velocities of about 2000 km/s for 
the clumps in the Vela supernova remnant (using their angu- 
lar distance to the pulsar, the age of the latter and the distance 
of the remnant), that are in reasonably good agreement with 
our Type II Model T310a. In fact, the Vela clumps show Mach 
cones and appear to move at present as sl ow as ~ 500 km/s 
(lAschenbachl2002tlAschenbach et alJ 19951) . indicating that su- 
personic drag plays an important role even in the supernova 
remnant phase. 
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